{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "This notebook gives suggests how to solve the problem of non-linear compressible flow using the automatic differentiation library included in PorePy. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model\n",
    "As an example, we will set up a non-linear problem for compressible flow. As usuall, we assume Darcy's law is valid:\n",
    "$$\n",
    "\\vec u = -\\mathcal K \\nabla p,\n",
    "$$\n",
    "where $\\vec u$ is the flux, $\\mathcal K$ the permeability tensor and $p$ the fulid pressure. Further, the conservation of mass gives\n",
    "$$\n",
    "\\frac{\\partial \\phi \\rho}{\\partial t} + \\nabla \\cdot \\rho \\vec u = q,\\quad \\text{in}\\ \\Omega \\\\\n",
    "u\\cdot n = 0,\\quad \\text{on}\\ \\partial \\Omega\n",
    "$$\n",
    "for porosity $\\phi$, fluid density $\\rho$, and source/sink term $q$.\n",
    "\n",
    "To solve this system of equation we need a constitutive law relating the fluid density to the pressure:\n",
    "$$\n",
    "\\rho = \\rho_r e^{c(p - p_r)},\n",
    "$$\n",
    "for reference density $\\rho_r$ and pressure $p_r$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import statements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.sparse as sps\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Porepy modules\n",
    "import porepy as pp\n",
    "import porepy.ad as ad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define constitutive laws and constants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set the porosity to 0.2 and let set the permeability to the default value (i.e. $\\mathcal K = 1$).\n",
    "We define the depenecy of $\\rho$ on $p$ as a function. Note that we have to use the exponent function ad.exp (and not np.exp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define data\n",
    "dt = 0.2                           # Time step\n",
    "phi = 0.2                          # Porosity \n",
    "c = 1e-1                           # Compressibility\n",
    "\n",
    "# Constitutive law\n",
    "def rho(p):\n",
    "    rho0 = 1\n",
    "    p_ref = 1\n",
    "    return rho0 * ad.exp(c * (p - p_ref))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discretization\n",
    "\n",
    "We use a finite-volume method to discretize the model equation. As a first step we create a partition of the domain into grid cells:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAFVBJREFUeJzt3Xu0ZGV55/Hvj26waWgFgsM0FwUj\nxmHIEBMUlJmEAaKoREiWQySjYi7Ta5KgaJLlqGNGVoyalUkcnRWTGYIoCQbEFpVcRlHUMReDCpqI\ntBHlDs3FC4KoI+Azf+zdbXFybn3sU+/p834/a9Xqqtq79vO8dU7Xr/bedepNVSFJ6tdurRuQJLVl\nEEhS5wwCSeqcQSBJnTMIJKlzBoEkdc4g0IqR5BtJHjfHshcl+Ztp9yT1wCDoTJKfS/Kp8UV3a5L/\nk+Tffh/bqySPn3HfhiRvTHJjkvuT3Jxkc5Jj5ttWVe1dVdcvsY89kpyT5Lqx5o1Jzk9y6FK2t1zG\nHi9s3cdySnLo+HuxtnUvWhyDoCNJfg14E/B64ADgMcAfAqcuYVuz/idP8gjgw8APA6cAjwT+FXAx\n8Mwd2dYO2gw8B/g54FHAUcBVwIk7YdvayQyJFaaqvHRwYXhx/AbwH+ZZ5ynAx4F7gK3AHwB7TCwv\n4FeB64AbgI+N990/bvtngV8aH7vXAv08bFsT9z1+vP4DwGXAvcAngNcCfzPHtk4CvgUcMk+9A8ft\nfRX4IvCfJpadA7wLuBC4D/gs8ATglcBdwC3A0yfW/yjwhrGve4H3AfuNy44Hbp1R+8axx5OB7wAP\njM/XP0z8bN46Pm+3Ab8NrJljHHsCFwBfA7YAL5+sN47z3cDd48/oJRPLHsHwRuD28fIm4BGTfY/b\nu2vs5TTgWcAXxuftVRPb2g14BfAl4CvAJRPPwc3jz/Ib4+WpwIuAvwX+x7j+68dt/vDENv8F8E3g\n0a3/v/R2cY+gH08F1gHvmWedh4CXAfuP658I/MqMdU4DjgGOqKofH+87qobDOu9keMH7QFXdv4ie\ntm9rlmVvAb4NbAR+YbzM5STgE1V1yzzrXMzwQncg8Fzg9UlOmFj+U8CfAvsCnwY+wPBidxDwW8D/\nnrG9F449bQQeBP7nPLUBqKr3M7wAvnN8vo4aF7193MbjgScBT2cI1Nm8BjgUeBzwk8Dzty1Ishvw\n58A/jH2fCLw0yTPGVf4rcCzwIwx7TE8BXj2x7X/J8DtyEPDfgD8et/9jwL8DfjPJYeO6L2b4+f0E\nw3P6NYafGcC234t9xnF+fLx9DHA9w97oaxl+Jtv7B84Arqiqu+cYu5ZL6yTyMp0L8B+BO3bwMS8F\n3jNxu4ATZqyz/V38ePtDwO9M3P4Rhj2Me4F/Wsy2gDUM75qfOLHs9cy9R/DHwMXzjOMQhpDbMHHf\nG4C3j9fPAT44seynGN7Jrhlvbxh722e8/dEZYzyC4Z3+GubZI5iodeHEsgOA/wfsOXHfGcBH5hjL\n9cAzJm7/0rZ6DC+0N89Y/5XA28brXwKeNbHsGcCN4/XjGfaqZo75mIn1rwJOG69vAU6cWLZx/Jmt\nZQiqAtZOLH/RLL0dw7D3kPH2p4DTW/9f6fHicbp+fAXYP8naqnpwthWSPAF4I3A0sJ7hP/VVM1ab\n7133tjobt92oqs8A+yQ5CThvkdt69Fh7cvlNC9R8wjzLDwS+WlX3zdje0RO375y4/i3gy1X10MRt\ngL0ZQo1ZetudYU9qRz12fOzWJNvu2425n5sDZyybvP5Y4MAk90zctwb464nHTj6PN433bfOVWcY8\n83nZe6LWe5J8d2L5QwzBNpeHjamqrkzyTeD4JFsZ3gRcNs/jtUw8NNSPjzO88zxtnnX+CPg8cHhV\nPRJ4FZAZ6yz0dbVXAE9PstciepprW3czHCo5ZOK+x8yznQ8BT0ly8BzLbwf2S7JhxvZuW0SPc5nZ\n2wPAlxnOl6zftiDJGoZg22bmmG9h+LnsX1X7jJdHVtW/nqPuVmBynJN93MJwvmWficuGqnrWuPx2\nhhfwyb5vn3eUc7sFeOaMWuuq6rZZxrjNbPdfwHB46AXA5qr69hL70ffBIOhEVX2d4bjvW5KclmR9\nkt2TPDPJ746rbWA4hPONJE8EfnkRm76T4Xj1Nn/C8GL1niRHJlmTZB0Pf/e9UK8PAZcC54x9HgGc\nOc/6HwI+ONb8sSRrx4+w/uckv1DDuYO/A96QZF2SfwP8IsPJ4aV6fpIjkqxnOIeweez7C8C6JM9O\nsjvDMfhHTDzuTuDQ8Xg+VbUVuBz4/SSPTLJbkh9M8hNz1L0EeGWSfZMcBJw1sewTwH1J/kuSPcfn\n/sgkTx6XXwS8Osmjk+zP8Puw1OfgfwGvS/JYgHGb2z59djfwXR7+ezGXC4GfZgiDP1liL/o+GQQd\nqarfB36N4cXpboZ3dWcB7x1X+Q2Gj1/ex3Dc/Z2L2Ow5wAVJ7kly+viO7t8D1wJ/yXhuAHgycPoO\ntHsWw2GIOxhOpr5tgfWfC/zV2PPXgWsYwudD4/IzGI5d385wwvw1Y4As1Z+Ofd3BcIL1JbA9cH+F\n4TDYbQx7CLdOPO5d479fSXL1eP2FwB4Mz9nXGD4Ku5HZ/da4vRvGsW1m2KPYFqCnMJyXuYFhD+U8\nhk8lwfBppE8B/8jwyairx/uW4s0Mh3EuT3If8PcMx/ypqm8CrwP+dvy9OHaujYwhfTXD3sJfz7We\nlte2kzSSFinJRxlO+M4859Gil18GnldVc+1BrHhJzgdur6pXL7iyloUni6VdSJKNDIdcPg4cDvw6\nw9977JLGv/z+GYaPzaoRDw1Ju5Y9GP6m4T6Gv+B+H8Nfh+9ykryW4RDef6+qG1r30zMPDUlS59wj\nkKTO7RLnCHZLmuy3hIU/NL/aajvm1V+3ZW3HPHVVVQu+4d8lgqDo7wfY4y9tb2P2ue6jduMxz/yD\n0Fl5aEiSOmcQSFLnDAJJ6pxBIEmdMwgkqXMGgSR1ziCQpM4ZBJLUOYNAkjpnEEhS5wwCSeqcQSBJ\nnTMIJKlzBoEkdc4gkKTOLVsQJDk/yV1Jrpm4b78kH0xy3fjvvstVX5K0OMu5R/B24OQZ970CuKKq\nDgeuGG9LkhpatiCoqo8BX51x96nABeP1C4DTlqu+JGlxpj1V5QFVtXW8fgdwwFwrJtkEbNp+e5kb\nm83aRnVb1nbMq79uy9qOeWVqNmdxVVWSOafyrKpzgXMB4uT1q75uy9q91W1Z2zFPv/ZiTPtTQ3cm\n2Qgw/nvXlOtLkmaYdhBcBpw5Xj8TeN+U60uSZkgt00GXJBcBxwP7A3cCrwHeC1wCPAa4CTi9qmae\nUJ5tWx4aWuV1W9burW7L2o65Qe2qBY8QLVsQ7EwGweqv27J2b3Vb1nbMDWovIgj8y2JJ6pxBIEmd\nMwgkqXMGgSR1ziCQpM4ZBJLUOYNAkjpnEEhS5wwCSeqcQSBJnTMIJKlzBoEkdc4gkKTOGQSS1Llm\nU1XuiNDfXKM9zq/a25h9rvuo7ZzFO0nR6feId1S3Ze3e6ras7ZinX3sxPDQkSZ0zCCSpcwaBJHXO\nIJCkzhkEktQ5g0CSOmcQSFLnDAJJ6pxBIEmdMwgkqXMGgSR1ziCQpM4ZBJLUOYNAkjpnEEhS55oE\nQZKXJflckmuSXJRkXYs+JEkNgiDJQcBLgKOr6khgDfC8afchSRq0OjS0FtgzyVpgPXB7oz4kqXtT\nn6qyqm5L8nvAzcC3gMur6vKZ6yXZBGzafnt6LW7n/Kp91O6tbsvajnllStV0Z9NMsi/wbuBngXuA\ndwGbq+rCeR4z5S7Huji/ag+1e6vbsrZjblC7asEcanFo6CTghqq6u6oeAC4FntagD0kSbYLgZuDY\nJOuTBDgR2NKgD0kSDYKgqq4ENgNXA58dezh32n1IkgZTP0ewFJ4jWP11W9burW7L2o65Qe0Veo5A\nkrSCGASS1DmDQJI6ZxBIUucMAknqnEEgSZ0zCCSpcwaBJHXOIJCkzhkEktQ5g0CSOmcQSFLnDAJJ\n6tzUp6pcitDfFHM9TqvX25h9rvuovStMVblLBEHR6dfHdlS3Ze3e6ras7ZinX3sxPDQkSZ0zCCSp\ncwaBJHXOIJCkzhkEktQ5g0CSOmcQSFLnDAJJ6pxBIEmdMwgkqXMGgSR1ziCQpM4ZBJLUOYNAkjpn\nEEhS55oEQZJ9kmxO8vkkW5I8tUUfkqR2E9O8GXh/VT03yR7A+kZ9SFL3UjXduXOSPAr4DPC4WmTx\nJFPucqyLsyn1ULu3ui1rO+YGtasWnKisxR7BYcDdwNuSHAVcBZxdVfdPrpRkE7Bp++2ptjhwftU+\navdWt2Vtx7wytdgjOBr4e+C4qroyyZuBe6vqN+d5jHsEq7xuy9q91W1Z2zE3qL2IPYIWJ4tvBW6t\nqivH25uBH23QhySJBkFQVXcAtyT5ofGuE4Frp92HJGnQ6lNDLwbeMX5i6Hrg5xv1IUndm/o5gqXw\nHMHqr9uydm91W9Z2zA1qr9BzBJKkFcQgkKTOGQSS1LkFgyDJi5PsO41mJEnTt5g9ggOATya5JMnJ\nSVb6H8lJknbAoj41NL74P53hY55HA5cAb62qLy1ve9vr+6mhVV63Ze3e6ras7Zgb1N5Znxoavxzu\njvHyILAvsDnJ734/TUqS2ltwjyDJ2cALgS8D5wHvraoHkuwGXFdVP7jsTbpHsOrrtqzdW92WtR1z\ng9o76dtH9wN+pqpumryzqr6b5JQl9idJWiH8y+L56uK7lx5q91a3ZW3H3KC2f1ksSVqIQSBJnTMI\nJKlzBoEkda7VfAQ7JPQ312iP86v2Nmaf6z5q7wpzFu8SQVB0era/o7ota/dWt2Vtxzz92ovhoSFJ\n6pxBIEmdMwgkqXMGgSR1ziCQpM4ZBJLUOYNAkjpnEEhS5wwCSeqcQSBJnTMIJKlzBoEkdc4gkKTO\nGQSS1DmDQJI61ywIkqxJ8ukkf9GqB0lS2z2Cs4EtDetLkmgUBEkOBp4NnNeiviTpe1pNVfkm4OXA\nhrlWSLIJ2LT99hSamsn5Vfuo3VvdlrUd88o09SBIcgpwV1VdleT4udarqnOBc8fHVI9zjTpm6662\n2o55+rUXo8WhoeOA5yS5EbgYOCHJhQ36kCQBqSbvtcfiwx7Bb1TVKQus5x7BKq/bsnZvdVvWdswN\nalctuGPg3xFIUuea7hEslnsEq79uy9q91W1Z2zE3qO0egSRpIQaBJHXOIJCkzhkEktQ5g0CSOmcQ\nSFLnDAJJ6pxBIEmdMwgkqXMGgSR1ziCQpM4ZBJLUOYNAkjrXaqrKHRL6m2Kux2n1ehuzz3UftZ2q\ncicpOv362I7qtqzdW92WtR3z9GsvhoeGJKlzBoEkdc4gkKTOGQSS1DmDQJI6ZxBIUucMAknqnEEg\nSZ0zCCSpcwaBJHXOIJCkzhkEktQ5g0CSOmcQSFLnDAJJ6tzUgyDJIUk+kuTaJJ9Lcva0e5AkfU+q\npjtlQpKNwMaqujrJBuAq4LSqunaex0y5y7EuTqLRQ+3e6ras7Zgb1K5acH6aqe8RVNXWqrp6vH4f\nsAU4aNp9SJIGTaeqTHIo8CTgylmWbQI2bb89ta6+x/lV+6jdW92WtR3zyjT1Q0PbCyd7A/8XeF1V\nXbrAuh4aWuV1W9burW7L2o65Qe2VeGgIIMnuwLuBdywUApKk5dXiU0MB3gpsqao3Tru+JOnhWuwR\nHAe8ADghyWfGy7Ma9CFJouE5gh3hOYLVX7dl7d7qtqztmBvUXqnnCCRJK4dBIEmdMwgkqXMGgSR1\nziCQpM4ZBJLUOYNAkjpnEEhS5wwCSeqcQSBJnTMIJKlzBoEkdc4gkKTOGQSS1LmmcxYvVuhvrtEe\n51ftbcw+133U3hXmLN4lgqDo9HvEO6rbsnZvdVvWdszTr70YHhqSpM4ZBJLUOYNAkjpnEEhS5wwC\nSeqcQSBJnTMIJKlzBoEkdc4gkKTOGQSS1DmDQJI6ZxBIUucMAknqnEEgSZ0zCCSpc02CIMnJSf4p\nyReTvKJFD5KkQaqmO2VCkjXAF4CfBG4FPgmcUVXXzvOYKXc51sVJNHqo3VvdlrUdc4PaVQvOT9Ni\nj+ApwBer6vqq+g5wMXBqgz4kSbSZqvIg4JaJ27cCx8xcKckmYNP228vf16xazjXqmK27Gms75qla\n1M7Iip2zuKrOBc4FSPKpqjq6cUtT5ZhXv97GC455pWpxaOg24JCJ2weP90mSGmgRBJ8EDk9yWJI9\ngOcBlzXoQ5JEg0NDVfVgkrOADwBrgPOr6nMLPOzc5e9sxXHMq19v4wXHvCJN/eOjkqSVxb8slqTO\nGQSS1LkVHQS9fRVFkkOSfCTJtUk+l+Ts1j1NS5I1ST6d5C9a9zINSfZJsjnJ55NsSfLU1j0ttyQv\nG3+vr0lyUZJ1rXva2ZKcn+SuJNdM3Ldfkg8muW78d9+WPc5mxQbB+FUUbwGeCRwBnJHkiLZdLbsH\ngV+vqiOAY4Ff7WDM25wNbGndxBS9GXh/VT0ROIpVPvYkBwEvAY6uqiMZPijyvLZdLYu3AyfPuO8V\nwBVVdThwxXh7RVmxQUCHX0VRVVur6urx+n0MLw4Hte1q+SU5GHg2cF7rXqYhyaOAHwfeClBV36mq\ne9p2NRVrgT2TrAXWA7c37menq6qPAV+dcfepwAXj9QuA06ba1CKs5CCY7asoVv2L4jZJDgWeBFzZ\ntpOpeBPwcuC7rRuZksOAu4G3jYfDzkuyV+umllNV3Qb8HnAzsBX4elVd3rarqTmgqraO1+8ADmjZ\nzGxWchB0K8newLuBl1bVva37WU5JTgHuqqqrWvcyRWuBHwX+qKqeBNzPCjxcsDONx8VPZQjBA4G9\nkjy/bVfTV8Pn9VfcZ/ZXchB0+VUUSXZnCIF3VNWlrfuZguOA5yS5keHw3wlJLmzb0rK7Fbi1qrbt\n7W1mCIbV7CTghqq6u6oeAC4Fnta4p2m5M8lGgPHfuxr388+s5CDo7qsokoThuPGWqnpj636moape\nWVUHV9WhDD/jD1fVqn6nWFV3ALck+aHxrhOBOefjWCVuBo5Nsn78PT+RVX6CfMJlwJnj9TOB9zXs\nZVYr+dtHl/JVFLu644AXAJ9N8pnxvldV1V817EnL48XAO8Y3OdcDP9+4n2VVVVcm2QxczfDpuE+z\nC3z1wo5KchFwPLB/kluB1wC/A1yS5BeBm4DT23U4O79iQpI6t5IPDUmSpsAgkKTOGQSS1DmDQJI6\nZxBIUucMAknqnEEgSZ0zCKQlSPLkJP+YZF2Svcbv2T+ydV/SUvgHZdISJfltYB2wJ8N3B72hcUvS\nkhgE0hKNXw/xSeDbwNOq6qHGLUlL4qEhael+ANgb2MCwZyDtktwjkJYoyWUMX519GLCxqs5q3JK0\nJCv220ellSzJC4EHqurPxvm1/y7JCVX14da9STvKPQJJ6pznCCSpcwaBJHXOIJCkzhkEktQ5g0CS\nOmcQSFLnDAJJ6tz/B8PPgxqiJvM7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1a3c4087b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Create grid\n",
    "g = pp.CartGrid([11,11])\n",
    "g.compute_geometry()\n",
    "pp.plot_grid(g, plot_2d=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, the model equation is integrated over each controll volume (i.e., each cell of the grid), and the divergence theorem is applied to the flux term:\n",
    "$$\n",
    "\\int_\\Omega \\phi \\frac{\\partial \\rho}{\\partial t} dV - \\int_{\\partial\\Omega}\\vec n\\cdot(\\rho\\vec u)dS - \\int_\\Omega q dV= 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The key-point of the finite-volume discretization is how the flux-term $\\vec u$ is approximated. We do not cover that in this tutorial(see e.g., I. Aavatsmark. An introduction to multipoint flux approximations for quadrilateral grids. Comput. Geosci., Vol. 6, No. 3, pp. 405–432, 2002. DOI: 10.1023/A:1021291114475).\n",
    "However, the main idea is that the fluid flux $\\vec u$ across a face is expressed as a linear combination of the cell-centered pressures $\\vec u = \\text{flux}\\ \\vec p$. Here, $\\text{flux}$ is the discretization matrix and $\\vec p$ is the vector of all cell-centered pressures.\n",
    "\n",
    "In porepy we can obtain the discretization matrix with, e.g, the two-point flux approximation:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize default data (i.e., unit parameters)\n",
    "data = pp.initialize_default_data(g, {}, 'flow')\n",
    "# Define flux discretization:\n",
    "flx_disc = pp.Tpfa('flow')\n",
    "# Discretize\n",
    "flx_disc.discretize(g, data)\n",
    "# The flux discretization can now be found in the dictionary as:\n",
    "flux = data[pp.DISCRETIZATION_MATRICES]['flow']['flux']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the negative sign in front of the surface-integral is included into the flux discretization matrix.\n",
    "\n",
    "The density is defined at the cell centers, but in the flux term we need to evaluate it at the faces. To do so, we simply take the average of the two neighbooring cells (note that other alternatives, such as upstream weighting, are commonly used).\n",
    "\n",
    "We also create discretized versions of the divergence operator div. The discrete divergence operator sums the fluxes in and out of each grid cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_faces_T = g.cell_faces.T\n",
    "def div(x):\n",
    "    \"\"\"\n",
    "    Discrete divergence\n",
    "    \"\"\"\n",
    "    return cell_faces_T * x\n",
    "\n",
    "def avg(x):\n",
    "    \"\"\"\n",
    "    Averageing. Note that this is not strictly correct for the boundary faces since\n",
    "    these only have 1 cell neighboor, but we have zero flux condition on these, so \n",
    "    this is not a problem.\n",
    "    \"\"\"\n",
    "    return 0.5 * np.abs(g.cell_faces) * x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Residual function\n",
    "To discretize the time-deriveative, we use backward Euler. Further, we assume that the densities are constant over each cell so we can take them out of the integral:\n",
    "$$\n",
    "\\int_\\Omega \\phi \\frac{\\rho^k - \\rho^{k-1}}{\\Delta t} dV =\\phi \\frac{\\rho^k - \\rho^{k-1}}{\\Delta t} \\int_\\Omega dV = \\phi \\frac{\\rho^k - \\rho^{k-1}}{\\Delta t}V,\n",
    "$$\n",
    "where $V$ is the volume of the cell. The same is also done for the source term.\n",
    "\n",
    "This gives us the residual\n",
    "$$\n",
    "\\phi \\frac{\\rho^k - \\rho^{k-1}}{\\Delta t} V + \\text{div}(\\text{avg}(\\rho^k)\\text{flux } p^k) - q^k V= 0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(p, p0):\n",
    "    # darcy:\n",
    "    u = flux * p\n",
    "\n",
    "    # Source:\n",
    "    src = np.zeros(g.num_cells)\n",
    "    src[60] = 1\n",
    "\n",
    "    # Define residual function\n",
    "    time = phi * (rho(p) - rho(p0)) / dt * g.cell_volumes\n",
    "    advection = div(avg(rho(p)) * u)\n",
    "    lhs = time + advection\n",
    "    rhs = src * g.cell_volumes\n",
    "\n",
    "    return lhs - rhs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialize Ad variable\n",
    "To initialize an AD variable create an Ad_array(...) with values equal the initial value and jacobian equal the identity matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set initial condition\n",
    "p0 = np.zeros(g.num_cells)\n",
    "p = ad.Ad_array(p0, sps.diags(np.ones(p0.shape)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time loop\n",
    "We are now ready to set up the time loop. We will set up a simple Newton iteration to find the zero of the residual function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solving time step:  1\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/rbe051/anaconda3/envs/porepy/lib/python3.6/site-packages/mpl_toolkits/mplot3d/axes3d.py:738: UserWarning: Attempting to set identical bottom==top results\n",
      "in singular transformations; automatically expanding.\n",
      "bottom=0.0, top=0.0\n",
      "  'bottom=%s, top=%s') % (bottom, top))\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVAAAADxCAYAAACd3+8mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl8VOW9+PHPMzPZCWsCCUlkkQgE\nVMBgwGqvbIq9iixRQahYRFukLV5tXa5L0R9XsdcNS60XrzsvidZegZcsWtHaWwVLEG+LuLAkkAWR\nACFkmWSW7++PIWMSSTJzMslMyPf9ep3XK+fMeZYzM/nOc5bneYyIoJRSKni2cFdAKaU6Kw2gSill\nkQZQpZSySAOoUkpZpAFUKaUs0gCqlFIWaQBVSp3xjDEvGGO+NcbsauZ1Y4x52hiz1xjzD2PMmEDy\n1QCqlOoKXgKmtvD6FUDmqeUW4A+BZKoBVCl1xhORvwLHWtjlauAV8dkG9DTGpLaWryPYegS5v1Kq\n6zJtSTzEGKkOcN9D8DngbLBplYisCqK4NKCowXrxqW2HWkoUbABVSqkOUQ38NMB9l4JTRLLbsTqn\npQFUKRWRDB0aoEqAjAbr6ae2tUivgUaYbt26sX///tO+9tJLL3HxxRd3cI2UCg8bEBfgEgLrgRtO\n3Y0fB5wQkRZP3+vr2OW89tprZGdn061bN1JTU7niiiv429/+Zjk/Ywx79+5ttO3kyZPcfvvtDBw4\nkISEBM466yxyc3P55JNPWsyrsrKSwYMHW6pHXV0dS5cuJTMzk4SEBAYOHMiCBQsoLCy0lF97Wbp0\nKfPmzQt3NdpVYWEhxhjcbne4q9JpGSAqwKXVvIxZA2wFhhpjio0xNxljfmaM+dmpXTYC+4G9wHPA\nrYHUscudwj/xxBMsX76cZ599lssvv5zo6Gg2b97MunXrgm7dud1uHI7vv4W1tbVMnDiRnj178vbb\nbzN8+HCcTiebNm1i06ZN5OTkBJxXMHJzcykuLua1115j9OjRVFVVsXr1arZs2cJNN93UprxV6IXi\nMz+ThfIUXkTmtPK6AIutZBzM0qmVl5dLQkKCvPHGG83u88knn8i4ceOkR48ekpKSIosXL5ba2lr/\n64CsXLlShgwZIgMHDpRLLrlEAImPj5eEhATJy8uT5557TlJSUqSysrLF+jTNq37bnj17RESkrKxM\nrrrqKklMTJSxY8fKfffdJz/4wQ9Om9ef//xniY2NlYMHDzZbXklJiVx11VXSq1cvOfvss2XVqlX+\n137zm99Ibm6uzJ07V7p16yYjR46Ur776Sh5++GFJTk6W9PR0eeedd/z7/8u//IvcfffdMnbsWElM\nTJRp06bJ0aNHRUTkgw8+kLS0tEZlDxgwQP785z/Lpk2bJCoqShwOhyQkJMh5550nIr7PZsGCBZKS\nkiL9+/eXe++9V9xu92mPo7q6Wm644Qbp2bOnDBs2TB599NFG5ZWUlMjMmTMlKSlJBg4cKCtWrPC/\n5nQ6ZcmSJZKamiqpqamyZMkScTqdjer96KOPSnJysqSkpMhbb70lGzZskMzMTOnVq5f8x3/8hz8v\nj8cjjzzyiAwePFh69+4t11xzjf89yMjIEEASEhIkISFBPv74Y3nxxRfloosukttuu0169+4t99xz\nj/Tq1Uv+8Y9/+PM8fPiwxMXFybffftvs59iJBBtfGi0DQJ4LcAHy21qelaVLBdBNmzaJ3W4Xl8vV\n7D75+fmydetWcblcUlBQIMOGDZMnn3zS/zogkydPlqNHj0p1dbV/W33QExG57rrrZP78+a3Wp7W8\nrrvuOrnmmmuksrJS/vnPf0r//v2bDaB33XWX/PCHP2yxvEsuuUQWLVokNTU1snPnTklKSpItW7aI\niC+AxsTEyObNm8XlcsmPf/xjGThwoCxbtkzq6upk1apV/iAv4gug/fv3l3/+859SWVkpM2fOlLlz\n54pIywG0vqz6fetNnz5dbrnlFqmsrJTDhw/L2LFj5dlnn23xWI8dOyZFRUVy7rnn+svzeDwyZswY\nefDBB6W2tlb27dsngwYNks2bN4uIyP333y85OTly+PBh+fbbb2X8+PFy3333+ettt9vlwQcf9B9z\nUlKSzJkzRyoqKmTXrl0SGxsr+/fvFxGRp556SnJycqSoqEicTqfccsstMnv2bBERKSgoEKDRd+3F\nF18Uu90uTz/9tLhcLqmurpZFixbJnXfe6d/nqaeekiuvvLLFz7ETaVNwGgjyYoCLBtAOsHr1aunX\nr19QaZ588kmZPn26fx3wB52G2xoG0EmTJsldd93lX9+5c6f06NFDEhMT5ZxzzgkoL7fbLQ6HQ774\n4gv/a/fcc0+zAXThwoVy3XXXNXscBw8eFJvNJhUVFf5td999tz/Q/+Y3v5HJkyf7X1u/fr0kJCT4\nW4EVFRUCyPHjx0XEF0AbHuPnn38uUVFR4na7gw6g33zzjURHR/t/REREXnvtNbn00ktPeywNA6KI\nyHPPPecvb9u2bZKRkdFo/4cfflhuvPFGEREZPHiwbNiwwf/a5s2bZcCAASLiC6CxsbHfO+Zt27b5\n9x8zZoy89dZbIiIybNgwee+99/yvlZaWisPh8P/4ni6ANq1bfX29Xq+IiFxwwQXy+uuvn/a4O6E2\nBadBIKsDXMIVQLvUBZg+ffpQVlbW4rWnr7/+mttvv538/Hyqq6txu91ccMEFjfbJyMg4bdqG5Rw6\n9N0NvFGjRlFeXs57773HwoULA8rryJEjuN3uRq8PGDCgxTK//vrrZl8vLS2ld+/eJCYmNsovPz/f\nv96vXz//33FxcSQlJWG32/3r4LvJ1bNnz+/VfcCAAbhcLsrKypqtQ3MOHDiAy+UiNfW7jh9er7fZ\n96a0tLTRaw3/PnDgAKWlpf46Ang8Hi655BJ/2obv44ABAygtLfWv9+nT53vH3PR9qays9Jc1Y8YM\nbLbv7sXa7XYOHz7c7LE2PaacnBzi4+P5y1/+QmpqKnv37mXatGnNpu9K6u/CR7IudRd+/PjxxMTE\nsHbt2mb3WbRoEcOGDWPPnj1UVFTw8MMP+5rqDRjTcgeLSZMm8e6771JVVdVqnZrLKzk5GYfDQVHR\nd50jDh482Gw+kydP5u9//zvFxcWnfb1///4cO3aMkydPNsovLS2t1To2p2ndoqKiSEpKIiEhgerq\n7/qQeDwejhw54l9veswZGRnExMRQVlZGeXk55eXlVFRU8Pnnn5+23NTU1EbH2bAeGRkZDBo0yJ9P\neXk5J0+eZOPGjYDvfThw4ECjevfv39/S8WdkZLBp06ZGZTmdTtLS0pr9XE+3ff78+axevZpXX32V\n3NxcYmNjLdXnTBPKu/DtpUsF0B49evDQQw+xePFi1q5dS3V1NS6Xi02bNnHnnXcCvsePunfvTrdu\n3fjyyy/5wx9aH1OgX79+jZ7dvOGGG0hNTWXGjBns2rULj8eD0+ls1Nprjd1uZ+bMmSxdupTq6mp2\n797Nyy+/3Oz+kydPZsqUKcyYMYMdO3bgdrs5efIkzz77LC+88AIZGRlcdNFF3HPPPTidTv7xj3/w\n/PPPt+lxotWrV7N7926qq6t54IEHyM3NxW63c8455+B0OtmwYQMul4tly5ZRW1vrT9evXz8KCwvx\ner2ALyBedtll3HHHHVRUVOD1etm3bx8ffvjhacu99tpreeSRRzh+/DglJSWsXLnS/9qFF15IYmIi\njz76KDU1NXg8Hnbt2sX27dsBmDNnDsuWLePIkSOUlZXx0EMPWX4Pfvazn3Hvvff6A/KRI0dYt24d\n4PsBtNlszT7T29C8efN46623WL16NTfccIOlupypHAEu4dKlAijAHXfcwRNPPMGyZctITk4mIyOD\nlStXMn36dAAee+wxXnvtNRITE7n55pu57rrrWs1z6dKlzJ8/n549e/LGG28QGxvLBx98QFZWFv/6\nr/9K9+7dGTp0KNu3b+eNN94IuK4rV66ksrKSlJQUbrzxRn7yk5+0uP+bb77Jj370I6677jp69OjB\nyJEjyc/PZ/LkyQCsWbOGwsJC+vfvz4wZM3jwwQf9r1nx4x//mBtvvJGUlBScTidPP/004PuheuaZ\nZ1i4cCFpaWkkJCSQnp7uT3fNNdcAvtPlMWN8o4a98sor1NXVkZWVRa9evcjNzW10GaShBx54gPT0\ndAYNGsTkyZPJzc0lJiYG8P3wvP3223z22WcMGjSIpKQkFi5cyIkTJwC47777yM7O5rzzzuPcc89l\nzJgx3HfffZaOf8mSJUybNo3LLruMxMRExo0b53/ONz4+nnvvvZcf/OAH9OzZk23btjWbT0ZGBmPG\njMEY47/UoDpHC9Q0PT1thQ4mogC49NJLmTdv3veu6YbDH/7wB/Ly8pptsXYGCxYsoH///ixbtizc\nVQmlNg0mMtQY+a8A950AO0T7wivVukOHDrF//37Gjx/Pnj17ePzxx/n5z38e7mpZVlhYyP/8z/+w\nc+fOcFclouhNJKXaQV1dHT/96U9JTExk4sSJXH311dx6a0A97yLO/fffz8iRI/n1r3/NoEGDwl2d\niKKn8EqprqxNp/AjjJE1Ae57vp7CK6XUd+pboJFMA6hSKiJ18HiglkR6/ZRSXZS2QJVSyiJD5N+F\n1wCqlIpIBogKNEKFadxqDaBKqYhkDAQ83rQGUKWU+o4xEGUPdy1apg/Sq5DYvn075513Hk6nk6qq\nKkaMGMGuXbvCXS3VidW3QANZwlZHfZBehcp9992H0+mkpqaG9PR07rnnnnBXSYVXmx6kz44ykt87\nwIK+Dc+D9BpAVcjU1dUxduxYYmNj+fjjj/0DE6suq20BNNpIfnKABZVqTyTVyR09epTKykpcLhdO\np5OEhIRwV0l1Zp3gSXptgaqQmTZtGrNnz6agoIBDhw41GuhYdUlta4HGGMlPb30/ALNfW6CqE3vl\nlVeIiori+uuvx+PxcNFFF/H+++8zceLEcFdNdVYGiPCrQNoCVUq1l7a1QOOM5Ac4wp/5QlugSin1\nHQPEhLsSLdMAqpSKTJ3gJlKEV08p1WVpAFVKqTaI8JtIGkCVUpFJW6BKKWWRBlCllLJI78IrpZRF\nnaAFqsPZtcHmzZsZOnQoQ4YMYfny5SHNu6ioiAkTJpCVlcWIESNYsWJFSPNvyOPxMHr0aK688sp2\nyb+8vJzc3FyGDRvG8OHD2bp1a8jLePLJJxkxYgQjR45kzpw5OJ3ONue5YMEC+vbty8iRI/3bjh07\nxpQpU8jMzGTKlCkcP368zeWoZtQH0ECWMNEAapHH42Hx4sVs2rSJ3bt3s2bNGnbv3h2y/B0OB48/\n/ji7d+9m27Zt/P73vw9p/g2tWLGC4cOHt0veAEuWLGHq1Kl8+eWX/N///V/IyyopKeHpp58mPz+f\nXbt24fF4yMvLa3O+N954I5s3b260bfny5UyaNIk9e/YwadKkkP9wqgbqu3IGsoSJBtAmXC4Xbreb\n1rq4/v3vf2fIkCEMHjyY6OhoZs+ezbp160JWj9TUVMaMGQNAYmIiw4cPp6SkJGT51ysuLmbDhg0s\nXLgw5HkDnDhxgr/+9a/cdNNNAERHR9OzZ8+Ql+N2u6mpqcHtdlNdXU3//v3bnOcPf/hDevduPCDl\nunXrmD9/PgDz589n7dq1bS5HNSPELVBjzFRjzFfGmL3GmLtP8/pZxpgPjDE7jTH/MMb8qLU8NYA2\n4fF4KCgoaDWIlpSUkJGR4V9PT09vlwAHUFhYyM6dO8nJyQl53rfddhu//e1vsdna56tQUFBAcnIy\nP/nJTxg9ejQLFy6kqqoqpGWkpaXxq1/9irPOOovU1FR69OjBZZddFtIy6h0+fJjU1FQAUlJSOHz4\ncLuUo/juJlIgS2tZGWMHfg9cAWQBc4wxWU12uw94Q0RGA7OBZ1rLVwPoaRQXF1NYWMi+fftabYm2\nt8rKSmbNmsVTTz1F9+7dQ5r322+/Td++fbngggtCmm9DbrebTz/9lEWLFrFz504SEhJCftp7/Phx\n1q1bR0FBAaWlpVRVVbF69eqQlnE6xhiMadN4GaoloW2BXgjsFZH9IlIH5AFXN9lHgPp/sh5AaWuZ\nagBtRnFxMUVFRezduxev1/u919PS0igqKmq0f1paWkjr4HK5mDVrFnPnzmXmzJkhzRvgo48+Yv36\n9QwcOJDZs2fz/vvvM2/evJCWkZ6eTnp6ur/1nJuby6effhrSMt577z0GDRpEcnIyUVFRzJw5k48/\n/jikZdTr168fhw4dAuDQoUP07du3XcpRBBtAk4wx+Q2WW5rklgYUNVgvPrWtoaXAPGNMMbAR+EVr\nVdQA2ori4mL27duH29143tSxY8eyZ88eCgoKqKurIy8vj2nTpoWsXBHhpptuYvjw4dx+++0hy7eh\nRx55xN/azsvLY+LEiSFvuaWkpJCRkcFXX30FwJYtW8jKanrm1DZnnXUW27Zto7q6GhFhy5Yt7XZT\nbNq0abz88ssAvPzyy1x9ddNGjAqpwANomYhkN1hWWShtDvCSiKQDPwJeNca0GCM1gAaguLiYgoKC\nRqf0DoeDXr16cfnllzN8+HCuvfZaRowYEbIyP/roI1599VWef/55Ro0axahRo9i4cWPI8m9o6tSp\n7ZJvvd/97ndceOGFnHfeeXz22Wf8+7//e0jzz8nJITc3lzFjxtC9e3e8Xi+33NK0ARK8OXPmMH78\neL766ivS09N5/vnnufvuu/nP//xPMjMzee+997j77u/di1ChEtq78CVARoP19FPbGroJeANARLYC\nsUBSi1XUAZUbczqdLT6neNVV06mqqggyVwfgbnWvzpOmI8uK5DR2wBNkGkhM7ElFRZd4frRtAyqn\nGMkP8IqSebzlAZWNMQ7ga2ASvsC5HbheRD5vsM8m4HUReckYMxzYAqRJC0Eywp/zjzy+4BlsS/BH\nwJ+DTDMF+CDINBOAvwWZ5mJgW5BpAMYBO4JMc4GFssZh7ZisvHdWPqPgzwpOnmz16RgFIe3KKSJu\nY8zPgXfw/fK9ICKfG2MeAvJFZD1wB/CcMebf8DUWb2wpeIIGUKVUpApxV04R2UiTXzwReaDB37uB\nHwSTpwZQpVRk6gR94SO8ekqpLksDqFJKtYGOSK+UUhZoC1QppSzSAZU7nyNHjuB0OomJidF+zkqF\nk7ZAO5+kpCRsNhtVVVXExsbicOhbpFRYaADtfIwxREdH43A4cDqduFwuYmNjtTWqVEer78oZwTSA\nNsNmsxEfH4/L5aKqqoqYmBiioqLwfaLB9iSx4+u1EmyaCRbSXGwhzbgg09SnC3YYPCtlWT0mK++d\nlc/ISq+iCI8KkUJboJ1fVFRUo9aor+9zJHfLtNK90spUIVn4uhYH4xwLZWVh7ZgitfsnBB+ouyiD\nbziPCKajMQXAGENcXBzR0dHhroo6Q7T3CFhnhE4wJ5K2QIOgN5RUqDSdrE6dRic4hdcWqFIRQqdR\nPg2d1lgpFQidRrmJTnAKrwFUqQih0yg3EeJpjdtDhF9hUKpr69LTKGtXTqVUqHS5aZT1JpJSqi26\n9DTKneAUXgOoUhGsS0+j3AkCqM7K2URrs3JOnDiF4GditDJ7YySn6ciyzrQ0AA5EXN/bOmfOHP7y\nl79QVlZGv379ePDBB5k+fTrXXnstBw8eZMCAAbzxxhvfu9EUwdo2K+cII/mvB1jQuS3PytleIvwK\nQyTyENndMguCTDMIOBpkGoA+QHWQaeItlNUHa8cUqd0/obl++mvWrDnt9i1btlgo4wzQCa6BRnj1\nlFJdlt6FV0opi7QFqpRSFmkAVUopizSAKqWUdRLhY09rAFVKRSSxQV2ED6isAbQJp9OJx+PBbo/w\nnz6lznBiwG0PtK+Pt13r0hwNoE243W5qa2vxer1ERUURFRWFzaYdtpTqaGIMnoAHMa9r17o0RwNo\nE926dSM+Ph4RweVyUVNTA+APph072Vuwk7Y58D1EHmyaPkGmqU8X3wFlWTkmqxPedcTkdaD/doHz\nRPiZoH6Szaif3jg6Ohqv14vL5aK6uhpfT6RI7lVkoXeQw0IPXbeB1CDTHTLBl+U2WOvxFKm9lyD4\nQN01CQZPhM9gquemAbDZbMTExJCQkBDuqqgzhE4q1zrB4MYe0BIuGkCVCoNAJpV78sknGTFiBCNH\njmTOnDk4nc4OqFnkEAx1xAS0hIsGUKUiUElJCU8//TT5+fns2rULj8dDXl5euKvVoepP4QNZwkUD\nqFIRyu12U1NTg9vtprq6mv79+4e7Sh0ulAHUGDPVGPOVMWavMebuZva51hiz2xjzuTHmtdby1JtI\nSkWgtLQ0fvWrX3HWWWcRFxfHZZddxmWXXRbuanWo+mugoWCMsQO/B6YAxcB2Y8x6EdndYJ9M4B7g\nByJy3BjT6vD/2gJVKgIdP36cdevWUVBQQGlpKVVVVaxevTrc1epQvlN4R0BLAC4E9orIfhGpA/KA\npsP73wz8XkSOA4jIt61lqgFUqQj03nvvMWjQIJKTk4mKimLmzJl8/PHH4a5Wh/LdRIoOaAGSjDH5\nDZZbmmSXBhQ1WC8+ta2hc4BzjDEfGWO2GWNafVRCT+GVikBnnXUW27Zto7q6mri4OLZs2UJ2dofP\nWBFWAsGcwpeFYEoPB5AJXAqkA381xpwrIuUtJVBKRZicnBxyc3MZM2YMDoeD0aNHc8stTRtVZzoT\n6Ol5IEqAjAbr6ae2NVQMfCK+CasKjDFf4wuo25utoU4q11jrk8pdBriDzNXK5GMOC+V0VJqOLKuj\n0oR/UrkzUJsmlRuWnSDP5w8PaN+LzY4WJ5UzxjiAr4FJ+ALnduB6Efm8wT5TgTkiMt8YkwTsBEaJ\nSLMTeWkLNGhuYFuQacYBu1vdq7EsLE3AZqWr5AQLv4sfGHgwyHS/sVDWB1a7f1qZvM7KZxTsdwF8\n3wcViFA94ykibmPMz4F38P3yvSAinxtjHgLyRWT9qdcuM8bsxvfL+OuWgidoAFVKRahQ94UXkY3A\nxibbHmjwtwC3n1oCogFUKRWRBENthE/LqQFUKRWROsNoTBpAlVIRSQOoUkq1QTiHqguEBlClVESS\n0D4H2i4iu3ZKqS5LT+GVUsoi31346HBXo0UaQJs4duwYlZWVGGNOuyilOoaewndCvXv3plu3bojI\naRffWxZsTxI7vl4rwbA4g6U7yCBvHL7ePsGyOXw9i9q9LAvHZOm9s/IZ2bHWq8jO1KlTA5rWo6vT\nU/hOqvkWpxtrszd+HWSac7A0G6WVmTKD7ZIJ8BvDUrkrqCRLzaPWun9aOSZLM3la+YyC/S4AXBBQ\n8CwvL2fhwoXs2rULYwwvvPAC48ePt1Be56TXQJVSli1ZsoSpU6fy5ptvUldXd2pa7a5DA6hSypIT\nJ07w17/+lZdeegmA6OhooqNbv6HywAMP0Lt3b2677TYA7r33Xvr27cuSJUvas7rtojN05dQR6ZWK\nQAUFBSQnJ/OTn/yE0aNHs3DhQqqqqlpNt2DBAl555RUAvF4veXl5zJs3r72r2y50Vk6llCVut5tP\nP/2URYsWsXPnThISEli+fHmr6QYOHEifPn3YuXMn7777LqNHj6ZPn2BvqEWOSA+gegqvVARKT08n\nPT2dnJwcAHJzcwMKoAALFy7kpZde4ptvvmHBggXtWc12FcpZOduLtkCVikApKSlkZGTw1VdfAbBl\nyxaysgJ7zGrGjBls3ryZ7du3c/nll7dnNdtViGflbBfaAlUqQv3ud79j7ty51NXVMXjwYF588cWA\n0kVHRzNhwgR69uyJ3R7ZLbjW6F14pZQlo0aNIj8/P+h0Xq+Xbdu28cc//rEdatVx6qc1jmQaQIPm\nwPdgfDDs+B66Drac+ODTHAqy146VHkWAzWHzPRjf7mVZOCZL752Vz8hO8N8FaM9/u927d3PllVcy\nY8YMMjMz262cjtAZroFqAA2aTioH4NVJ5YjESeWysrLYv39/u+XfkbQvvFJKtYFeA1VKKQu0K6dS\nSlmk10CVUsoi3134yO4LrwFUKRWR9BReKaXaQAOoUkpZoNdAlVLKos7wHKgOJtIKj8dDXV0dNTU1\nAY3HqFQgpk6dGu4qRLz6rpyBLOES2eE9DKqqqnA6nXg8HkQEu92O3W4nOjoam82GTip3ik4qh/VJ\n5RwBTyjn8XjIzs4mLS2Nt99+20JZnZeewndCNpsNh8NBTExMC5PK/S3IXC/G2kR0BUGmGUTQk6lJ\nfPBdJcEX1KxM9mapW6aVCeKsvHdWPqNgvwvg+z4EZsWKFQwfPpyKigoL5XR+egrfycTFxeFwOHQO\neBV2xcXFbNiwgYULF4a7KmHRGab0iOzwrlQXdtttt/Hb3/6WkydPhrsqYdEZngPVFqhSEejtt9+m\nb9++XHCBleHyzhxu7AEt4aItUKUi0EcffcT69evZuHEjTqeTiooK5s2bx+rVq8NdtQ7jxRbxXTm1\nBapUBHrkkUcoLi6msLCQvLw8Jk6c2KWCZ71QXgM1xkw1xnxljNlrjLm7hf1mGWPEGJPdWp7aAlVK\nRaRQXgM1xtiB3wNTgGJguzFmvYjsbrJfIrAE+CSQfLUFqlSEu/TSS7vcM6AAQkivgV4I7BWR/SJS\nB+QBV59mv/8HPAo4A8lUA6hSKkIFNa1xkjEmv8FyS5PM0oCiBuvFp7Z9V5oxY4AMEdkQaA31FF4p\nFZGCPIUvE5FWr1k2xxhjA54AbgwmnQbQoDkIpieJj5XZGx34escEm8bCTJ5Bd5U8lc7KbJmWumVa\nOKag3zsrn5Gd4L8LoP92gREMtaHr514CZDRYTz+1rV4iMBL4y6lONCnAemPMNBFpdm5p/SSD5gY+\nCDLNBCK7+2ewM1iCr6+5lS6WVmbLjNRumRcT/HcBfN8H1ZoQj8a0Hcg0xgzCFzhnA9f7yxI5ASTV\nrxtj/gL8qqXgCRpAlVIRLFR34UXEbYz5OfAOvlOHF0Tkc2PMQ0C+iKy3kq8GUKVURAp1V04R2Qhs\nbLLtgWb2vTSQPDWAKqUikmDweCO7L7wGUKVURBKvodYZ2V05NYAqpSKSiMHj1haoUkoFT9AAeibx\neDzhroJSXYaIwe2K7ACqXTkD4Ha7qaqqora2NtxVUWeI1iaVKyoqYsKECWRlZTFixAhWrFjRQTWL\nJAavxxHQEi7aAm2B2+2mtradknpsAAAQ5UlEQVQWYwyxsbHY7ZH9a6g6j9YmlXM4HDz++OOMGTOG\nkydPcsEFFzBlyhSysoKd+K4TE0BP4TsXEfEHTpvNRlxc3KnZOOvZCb4niZUufx3Z/TPYGSzr01np\nYmlhptGI7ZZp5bsAgfzbpaamkpqaCkBiYiLDhw+npKSkawVQrwFnZIeoyK5dGJSVleFyuU4TOOt5\ngD8HmesUIrv75+5W9/q+LODrINOcY6GsLCK3W+YEgv8ugO/7ELjCwkJ27txJTk6OhbI6OXe4K9Ay\nDaBNJCcnExcXF+5qKAVAZWUls2bN4qmnnqJ79+7hrk7H8g0IGtE0gCoVoVwuF7NmzWLu3LnMnDkz\n3NXpeBpAlVJWiAg33XQTw4cP5/bbbw93dcJDAFe4K9EyfYxJqQj00Ucf8eqrr/L+++8zatQoRo0a\nxcaNG1tPeCYRoDbAJUy0BapUBLr44osRkXBXI7z0FF4ppSzSAKqUUhZpAFVKKYs0gJ6J7AT7IHRk\n916y43tYPVh2fA/Gt3dZkdyryMp3AfTfLggaQM80HprMChCAHxHZvZe2BZkGYBzWeggFW9Y4IrdX\n0RSC/y6A7/ugWuUFnOGuRMs0gCqlIpOewiullEUaQJVSyiINoEop1QYaQJVSygJtgSqllEVeoCbc\nlWiZDiYSIBGhuro63NVQZ4jW5kQC37QfQ4cOZciQISxfvrwDahVhBN9Tg4EsYaIBNAAej4eqqiqi\noqLCXRV1hmhtTiSPx8PixYvZtGkTu3fvZs2aNezebWXmgE7OHeASJhpAW+FyuaipqSE+Pl4DqAqZ\n1lqgf//73xkyZAiDBw8mOjqa2bNns27dug6qXYSovwYawQFUr4G2wOl04vV6SUhIwBiDiBAfn0h1\ndbA9STqq+6eD4LswOvD19gmWA2uT3gVblpVj6qhumXas9SpyUFxczNSpU5ttiZaUlJCRkeFfT09P\n55NPPrFQViemN5E6J6/XS01NDQ6Hg/h438yTIoLNZuP111/DbrcTExPTrnWoL7+9W701NTVER0e3\n+5TNNTU1xMTENDNRX+hUV1cTGxvbruXUXw+3OtX1r3/9aw4cOIDT6WwxiHZ52pWz8zlx4oT/n8Ph\n8L09NpuNiooKvF6vvyXqdrffT2P9QLper5e6urp2K6e+DK/X265lgO+Yamo65pZqVVVVuwdq8AVr\nK+UsXbqUuXPnUlRURFFRET169GD8+PGNAmlaWhpFRUX+9eLiYtLS0kJS705FW6CdS0xMDPHx8dhs\nNkQEu93uD2Lx8fHt3lLzeDzU1tYSFxeHMaZdy3K73bjdbmJjY9u1HIC6ujqMMR1yHbm6upqYmJh2\n/6zqz1SstETXrl3LXXfdRVlZGU6nk2+//bZRa3Ts2LHs2bOHgoIC0tLSyMvL47XXXmuPw4hcneAU\nXm8iNVF/+lcfPJ1OJ06nk7i4uHb/hxQRnE4nsbGx7R48wRdAO+rGWH3LvSPExMRQW9v+E+XYbDbi\n4uJwOp14PME/S/Poo4+SnJxMfHw8RUVFlJSU+G8uORwOVq5cyeWXX87w4cO59tprGTFiRKgPIbLV\nTyoXyBImJsgv9Rk/SUt1dTVbt27Fbrfjcrmora3tsIDmcrkwxvgvHbS32tradr+WW8/r9eLxeDos\nYNfW1hIdHd0hn1v9pRarLfm7776b48ePU15eTmpqKmlpaWfKddE2vfkmKVuYlh/Yzi+aHSKS3Zby\nrNAA2oSI8Le//c1Si0Ipq+68807Ky8uJiYkhJiaGpKSkMyGIti2A9skW/jXAAPpq6wHUGDMVWIHv\n8Yn/FpHlTV6/HViI78LBEWCBiBxoMU8NoEqpdtK2ANo7W5gUYAB9s+UAaoyxA1/je1atGNgOzBGR\n3Q32mQB8IiLVxphFwKUicl1Lxeo1UKVUZAptV84Lgb0isl9E6oA84OpGxYl8ICL1/bW3AemtZap3\n4ZVSkSm4u/BJxpiGzdVVIrKqwXoaUNRgvRjIaSG/m4BNrRWqAVQpFZmCC6BlobqJZIyZB2QD/9La\nvnoKr1Q7O3bsGFOmTCEzM5MpU6Zw/Pjx0+738ssvk5mZSWZmJi+//LJ/+44dOzj33HMZMmQIv/zl\nL/2Pg/3xj39kxIgR2Gw28vMbXyt85JFHGDJkCEOHDuWdd97xb+9UIzyF9jGmEiCjwXr6qW2NGGMm\nA/cC00Sk1Wfh9CaSUu3szjvvJDY2lq1bt7Jz5066devGzp076dWrl3+fY8eOkZ2dzR133METTzzB\nwYMHWbFiBbfeeisXXnght956K48//jj79+9n8uTJrF27li+//JKKigqmTp1KYmIiQ4cO5Y033uDQ\noUNMmjSJvn374nK52Lt3Lx6Ph2+++YZx48ZRW1tLz5492bdvH2effTa7du1qr0Nv202k7tlCdoA3\nkT5o9SaSA99NpEn4Aud24HoR+bzBPqOBN4GpIrInkGK1BapUO1u3bh1lZWWMHz+ezMxMDhw4QN++\nfbn00kv9rdF33nmHSy65hMcff5zbb78dh8PB4sWLSUpKoqioiGeeeYY77riDbt26sX79erp168Yz\nzzzDm2++Sffu3UlOTubDDz+kb9++XHnlldx8882sWLGCwsJC3G43IsLIkSPp168fDoeD+++/n7q6\nOr788kv69+9PdnaHP0LZuvq+8IEsrRARN/Bz4B3gC+ANEfncGPOQMWbaqd3+E+gG/NEY85kxZn1r\n+WoLVKl21rNnT/r160dmZibvv/++f3jE7Oxsxo0bx6OPPspjjz3G1q1b2bFjByUlJYgI55xzDqWl\npTidThISEqiqqsLj8ZCcnIzdbqdfv34UFhZSXl6O3W7noosuYteuXURFRTFs2DAOHTpEQUEBdrud\n66+/nr59+7Jq1Sqqq6upq6tDRJg+fTp79uxhzZo1ZGVlhfrQ29YCTcgWsgJsgeaH50F6bYEqFQKT\nJ09m5MiR31vqx/D85ptvePfdd/1DI1ZXV/PZZ5+xdu1afx5ffPEFNpuN3r17Y4xhz5499OvXD7vd\nTlVVFfHx8YgI33zzDVFRUXTv3p2jR49is9nweDx89NFHnDhxgqqqKrZu3YrD4cBut1NdXc0HH3zA\noEGDOHHihH8MBI/Hw9q1a0lOTo7MsUZ1RHqluob33nuPlJSU722/9957SUhIwO124/V6OXr0KMnJ\nyQBUVFRQWFgI+EZfKi0tRUSIiooiMTERt9vN/v37qa2txeFw4HQ6/cG1oKCA/fv3A74uwA6Hg+7d\nu+NyuaisrMTpdDJnzhz/QDj79u3jl7/8JQAzZ84kPT2d6OhoPB4PH374IX/605864F2yIMIHVNYA\nqlSI1AfRwsJC9u7d618OHz5MTU2N/+650+mke/fugG+Uqn/7t3/j8ssvp7KykrKyMmw2G5WVldjt\ndn9LsaamhsTERI4dO0ZSUhLGGA4cOODv619XV8fll18O4B+e8OabbyYpKQmbzYbNZqOmpgaXy8XH\nH39MTU2Nfyi+6OhoduzYwYIFCzr6LWtZJxiRXq+BKtXOjh49yrBhw/zBMS0tjSNHjuByufB4PNjt\ndn7xi1+Ql5fHN998gzGG1NRUSktLsdvtjcZliI6O9s+ScOLEiUbl2Gw2oqOjcTp9d1X69+9PaWkp\nM2bM4K233iIqKgq3240xBq/XS3x8PNXV1QwYMICDBw+SnJzMpk2bGDNmTKgOvW3XQGOyhbQAr4EW\n6DVQpc5Iffr04dVXX/WPRHX06FFiY2MREZKSkujVqxd5eXkMHjyYPn36ICKUlpbicDj8rcn6lqbL\n5SIqKorKykri4uL8I3cZY7DZbPTr149u3boBMGrUKNLT01m7dq0/vc1mIzExEcA/etTBgweJj49n\n7ty5LFq0qOPemNZ0ghaoBlClOsDkyZPp06cPMTExOJ1OysvLiYqK4qKLLmLw4ME4nU6Ki4uJi4vz\nn6LXP35U/7yo3W4nKiqKmpoavF4vycnJ/lZnXFwcbreb4uJioqKisNvt7N+/nyNHjiAi/ilpjDFc\nddVVREdH+x+hEhE8Hg85OTmUl5dz6NChcL5VjWkAVUo5HA5WrVpFXV0ddrudgQMHEh8f77+jnp2d\nzRVXXEFSUhLHjx8nOjqa5ORkrrjiCn+rsf5UPisri6SkJIYOHcr5559PXV0dtbW1jB49mhEjRjB9\n+nTi4uIoKfF1tJkxYwY9evTg/PPP59xzz+X48ePU1dUxYMAARo0aRd++fenduzf33HMP6enp/nRh\n1wkGVNYAqlQHueqqq3jrrbfwer2UlJQwcOBAnE6nf+psh8PB1q1bOffcc3G73VRUVJCTk8Oll17q\nf2zJ7Xazb98+pk+fzt69e+nZsydjx45FRPjiiy+YOXMmW7duJSYmhsrKSvr168f+/fs5efIkn376\nKRMnTmTq1KkYYzh69CgFBQXMnz+ftLQ0vF4vLlcYo1FT+hiTUqqhhkG0PuDVB9GKigp/l8/6IPr6\n669zySWX+IPorbfeitvt5k9/+hP333+/P4guX74ct9vN888/77+j/89//pOzzz4bj8fjvzG1dOlS\nZs2aRWZmJrNmzcLr9bJmzRoWL16M1+vl8OHDkTN5XSe4Buq/PhLgopQKgfXr14vD4ZCMjAxZunSp\nnHfeefLTn/5U1q1bJyIiNTU1MmrUKOnevbuMHTtWnnrqKbnmmmtEROQXv/iFREdHy5AhQ+SFF16Q\nQYMGidvtlg0bNkhmZqZkZGRI3759G5WXnp4u48ePl/PPP19SUlKke/fuMmTIEMnIyJDExES58MIL\n5b/+679k7NixoTzMYONLowUuEBwS2AL5bS3PUh2DTKCUCpH6gDd48GBZtmyZiIjcf//9jYJobm6u\nnH322TJ27FjZt2+fP+2yZctk8ODBcs4558jGjRv922fPni0pKSnicDgkLS1N/vu//1tERMrKymTi\nxIkyZMgQmThxouTk5Mj69evF6/XKrbfeKoMHD5aRI0fK9u3bQ3mIbQ+gRgJbwhRA9TlQpbqYW2+9\nleTkZB588MH2Lqptz4GabIEAnwMlPM+B6oDKSnUhL730EgcOHGDlypXhrsoZQQOoUl3Ejh07eOyx\nx/jf//1ffzdO1Tb6LirVRaxcuZJjx44xYcIERo0axcKFC8NdpU5Pr4EqpdpLG6+BjhH4KMC94/Ua\nqFJKfae+K1Lk0gCqlIpQwU3LGQ4aQJVSEUpboEopZZEGUKWUskiAmnBXokUaQJVSEUqvgSqllEV6\nCq+UUhZpC1QppSzSFqhSSlmkLVCllLLIi96FV0opS/QUXiml2kBP4ZVSygJtgSqllEUaQJVSyiK9\nC6+UUhbpXXillLJIT+GVUsqiyD+F10nllFIRqr4FGsjSOmPMVGPMV8aYvcaYu0/zeowx5vVTr39i\njBnYWp4aQJVSEaq+BRrI0jJjjB34PXAFkAXMMcZkNdntJuC4iAwBngQebS1fDaBKqQhVfxMpkKVV\nFwJ7RWS/iNQBecDVTfa5Gnj51N9vApOMMS3OLBrsNdA2TVOqlFKBO/QOLE0KcOdYY0x+g/VVIrKq\nwXoaUNRgvRjIaZKHfx8RcRtjTgB9gLLmCtWbSEqpiCQiU8Ndh9boKbxSqisoATIarKef2nbafYwx\nDqAHcLSlTDWAKqW6gu1ApjFmkDEmGpgNrG+yz3pg/qm/c4H3RURaylRP4ZVSZ7xT1zR/DrwD2IEX\nRORzY8xDQL6IrAeeB141xuwFjuELsi0yrQRYpZRSzdBTeKWUskgDqFJKWaQBVCmlLNIAqpRSFmkA\nVUopizSAKqWURRpAlVLKov8PDZukEliYvWkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f19ccef41d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solving time step:  2\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVAAAADxCAYAAACd3+8mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl8lOW58PHfPUtWwppAQhJZJLIr\nYDBia4+yCPYoCkQFoWIRbZG2WG0tHpeiL6diT12w1PbQ484r0doj8MpiK9r2VMESxFMRFwIJZkEk\nQAhZJpnlev8YMiaRJDNPJpkJub6fz/P5MDPPvTzPDFfuZ7me24gISimlQmeLdAeUUqqr0gCqlFIW\naQBVSimLNIAqpZRFGkCVUsoiDaBKKWWRBlCl1FnPGPOMMeZLY8zeFj43xpgnjTEFxph/GmMmBFOv\nBlClVHfwHDCjlc+vBLJOL7cBvw2mUg2gSqmznoj8DTjeyirXAC+I306gtzEmra16HaH2I8T1lVLd\nl2lP4WHGSE2Q6x6GjwBXo7fWisjaEJpLB4obvS45/d7h1gqFGkCVUqpT1ADfC3LdFeASkewO7M4Z\naQBVSkUlQ6cGqFIgs9HrjNPvtUrPgUaZHj16cPDgwTN+9txzz/HNb36zk3ukVGTYgPgglzDYBNx0\n+mr8xcBJEWn18L2hj93OSy+9RHZ2Nj169CAtLY0rr7ySv//975brM8ZQUFDQ5L1Tp05x5513Mnjw\nYBITEznnnHPIzc3lvffea7Wuqqoqhg4daqkf9fX1rFixgqysLBITExk8eDCLFi2iqKjIUn0dZcWK\nFSxYsCDS3ehQRUVFGGPweDyR7kqXZQBnkEubdRmzHtgBDDfGlBhjbjHGfN8Y8/3Tq2wBDgIFwO+B\n24PpY7c7hH/sscdYtWoVv/vd75g+fToxMTFs27aNjRs3hjy683g8OBxf34V1dXVMnjyZ3r178/rr\nrzNy5EhcLhdbt25l69at5OTkBF1XKHJzcykpKeGll15i/PjxVFdXs27dOrZv384tt9zSrrpV+IXj\nOz+bhfMQXkTmtfG5AEutVBzK0qVVVFRIYmKivPLKKy2u895778nFF18svXr1ktTUVFm6dKnU1dUF\nPgdkzZo1MmzYMBk8eLBceumlAkhCQoIkJiZKXl6e/P73v5fU1FSpqqpqtT/N62p4b//+/SIiUl5e\nLldffbUkJSXJxIkT5b777pNvfOMbZ6zrz3/+s8TFxcnnn3/eYnulpaVy9dVXS58+feTcc8+VtWvX\nBj77+c9/Lrm5uTJ//nzp0aOHjBkzRj799FP5xS9+ISkpKZKRkSFvvPFGYP1/+Zd/keXLl8vEiRMl\nKSlJZs6cKceOHRMRkbffflvS09ObtD1o0CD585//LFu3bhWn0ykOh0MSExPl/PPPFxH/d7No0SJJ\nTU2VgQMHyr333isej+eM21FTUyM33XST9O7dW0aMGCGPPPJIk/ZKS0tl9uzZkpycLIMHD5bVq1cH\nPnO5XLJs2TJJS0uTtLQ0WbZsmbhcrib9fuSRRyQlJUVSU1Pltddek82bN0tWVpb06dNH/v3f/z1Q\nl9frlYcffliGDh0qffv2leuuuy6wDzIzMwWQxMRESUxMlHfffVeeffZZueSSS+SOO+6Qvn37yj33\n3CN9+vSRf/7zn4E6jxw5IvHx8fLll1+2+D12IaHGlybLIJDfB7kA+e1tz8rSrQLo1q1bxW63i9vt\nbnGd/Px82bFjh7jdbiksLJQRI0bI448/HvgckKlTp8qxY8ekpqYm8F5D0BMRueGGG2ThwoVt9qet\num644Qa57rrrpKqqSj788EMZOHBgiwH0Zz/7mXzrW99qtb1LL71UlixZIrW1tbJnzx5JTk6W7du3\ni4g/gMbGxsq2bdvE7XbLd77zHRk8eLCsXLlS6uvrZe3atYEgL+IPoAMHDpQPP/xQqqqqZPbs2TJ/\n/nwRaT2ANrTVsG6Da6+9Vm677TapqqqSI0eOyMSJE+V3v/tdq9t6/PhxKS4ulrFjxwba83q9MmHC\nBHnwwQelrq5ODhw4IEOGDJFt27aJiMj9998vOTk5cuTIEfnyyy9l0qRJct999wX6bbfb5cEHHwxs\nc3JyssybN08qKytl7969EhcXJwcPHhQRkSeeeEJycnKkuLhYXC6X3HbbbTJ37lwRESksLBSgyW/t\n2WefFbvdLk8++aS43W6pqamRJUuWyN133x1Y54knnpCrrrqq1e+xC2lXcBoM8myQiwbQTrBu3ToZ\nMGBASGUef/xxufbaawOvgUDQafxe4wA6ZcoU+dnPfhZ4vWfPHunVq5ckJSXJeeedF1RdHo9HHA6H\nfPzxx4HP7rnnnhYD6OLFi+WGG25ocTs+//xzsdlsUllZGXhv+fLlgUD/85//XKZOnRr4bNOmTZKY\nmBgYBVZWVgogJ06cEBF/AG28jR999JE4nU7xeDwhB9AvvvhCYmJiAn9EREReeuklueyyy864LY0D\noojI73//+0B7O3fulMzMzCbr/+IXv5Cbb75ZRESGDh0qmzdvDny2bds2GTRokIj4A2hcXNzXtnnn\nzp2B9SdMmCCvvfaaiIiMGDFC3nzzzcBnZWVl4nA4An98zxRAm/etob8+n09ERC688EJ5+eWXz7jd\nXVC7gtMQkHVBLpEKoN3qBEy/fv0oLy9v9dzTZ599xp133kl+fj41NTV4PB4uvPDCJutkZmaesWzj\ndg4f/uoC3rhx46ioqODNN99k8eLFQdV19OhRPB5Pk88HDRrUapufffZZi5+XlZXRt29fkpKSmtSX\nn58feD1gwIDAv+Pj40lOTsZutwdeg/8iV+/evb/W90GDBuF2uykvL2+xDy05dOgQbrebtLSvEj98\nPl+L+6asrKzJZ43/fejQIcrKygJ9BPB6vVx66aWBso3346BBgygrKwu87tev39e2ufl+qaqqCrQ1\na9YsbLavrsXa7XaOHDnS4rY236acnBwSEhL4y1/+QlpaGgUFBcycObPF8t1Jw1X4aNatrsJPmjSJ\n2NhYNmzY0OI6S5YsYcSIEezfv5/Kykp+8Ytf+IfqjRjTeoLFlClT+NOf/kR1dXWbfWqprpSUFBwO\nB8XFXyVHfP755y3WM3XqVP7xj39QUlJyxs8HDhzI8ePHOXXqVJP60tPT2+xjS5r3zel0kpycTGJi\nIjU1X+WQeL1ejh49GnjdfJszMzOJjY2lvLyciooKKioqqKys5KOPPjpju2lpaU22s3E/MjMzGTJk\nSKCeiooKTp06xZYtWwD/fjh06FCTfg8cONDS9mdmZrJ169YmbblcLtLT01v8Xs/0/sKFC1m3bh0v\nvvgiubm5xMXFWerP2SacV+E7SrcKoL169eKhhx5i6dKlbNiwgZqaGtxuN1u3buXuu+8G/Lcf9ezZ\nkx49evDJJ5/w29+2/UyBAQMGNLl386abbiItLY1Zs2axd+9evF4vLperyWivLXa7ndmzZ7NixQpq\namrYt28fzz//fIvrT506lWnTpjFr1ix2796Nx+Ph1KlT/O53v+OZZ54hMzOTSy65hHvuuQeXy8U/\n//lPnn766XbdTrRu3Tr27dtHTU0NDzzwALm5udjtds477zxcLhebN2/G7XazcuVK6urqAuUGDBhA\nUVERPp8P8AfEK664grvuuovKykp8Ph8HDhzgr3/96xnbvf7663n44Yc5ceIEpaWlrFmzJvDZRRdd\nRFJSEo888gi1tbV4vV727t3Lrl27AJg3bx4rV67k6NGjlJeX89BDD1neB9///ve59957AwH56NGj\nbNy4EfD/AbTZbC3e09vYggULeO2111i3bh033XSTpb6crRxBLpHSrQIowF133cVjjz3GypUrSUlJ\nITMzkzVr1nDttdcC8Ktf/YqXXnqJpKQkbr31Vm644YY261yxYgULFy6kd+/evPLKK8TFxfH2228z\natQo/vVf/5WePXsyfPhwdu3axSuvvBJ0X9esWUNVVRWpqancfPPNfPe73211/VdffZVvf/vb3HDD\nDfTq1YsxY8aQn5/P1KlTAVi/fj1FRUUMHDiQWbNm8eCDDwY+s+I73/kON998M6mpqbhcLp588knA\n/4fqqaeeYvHixaSnp5OYmEhGRkag3HXXXQf4D5cnTPA/NeyFF16gvr6eUaNG0adPH3Jzc5ucBmns\ngQceICMjgyFDhjB16lRyc3OJjY0F/H94Xn/9dT744AOGDBlCcnIyixcv5uTJkwDcd999ZGdnc/75\n5zN27FgmTJjAfffdZ2n7ly1bxsyZM7niiitISkri4osvDtznm5CQwL333ss3vvENevfuzc6dO1us\nJzMzkwkTJmCMCZxqUF1jBGqaH562QR8mogC47LLLWLBgwdfO6UbCb3/7W/Ly8locsXYFixYtYuDA\ngaxcuTLSXQmndj1MZLgx8p9Brns57BbNhVeqbYcPH+bgwYNMmjSJ/fv38+ijj/KDH/wg0t2yrKio\niP/+7/9mz549ke5KVNGLSEp1gPr6er73ve+RlJTE5MmTueaaa7j99qAy76LO/fffz5gxY/jpT3/K\nkCFDIt2dqKKH8Eqp7qxdh/CjjZH1Qa57gR7CK6XUVxpGoNFMA6hSKip18vNALYn2/imluikdgSql\nlEWG6L8KrwFUKRWVDOAMNkJF6LnVGkCVUlHJGAj6edMaQJVS6ivGgNMe6V60Tm+kV2Gxa9cuzj//\nfFwuF9XV1YwePZq9e/dGuluqC2sYgQazRKyPeiO9Cpf77rsPl8tFbW0tGRkZ3HPPPZHukoqsdt1I\nn+00kt83yIa+jMyN9BpAVdjU19czceJE4uLiePfddwMPJlbdVvsCaIyR/JQgGyrTTCTVxR07doyq\nqircbjcul4vExMRId0l1ZV3gTnodgaqwmTlzJnPnzqWwsJDDhw83edCx6pbaNwKNNZKf0fZ6AOag\njkBVF/bCCy/gdDq58cYb8Xq9XHLJJbz11ltMnjw50l1TXZUBovwskI5AlVIdpX0j0Hgj+UE+4c98\nrCNQpZT6igFiI92J1mkAVUpFpy5wESnKu6eU6rY0gCqlVDtE+UUkDaBKqeikI1CllLJIA6hSSlmk\nV+GVUsqiLjAC1cfZtcO2bdsYPnw4w4YNY9WqVWGtu7i4mMsvv5xRo0YxevRoVq9eHdb6G/N6vYwf\nP56rrrqqQ+qvqKggNzeXESNGMHLkSHbs2BH2Nh5//HFGjx7NmDFjmDdvHi6Xq911Llq0iP79+zNm\nzJjAe8ePH2fatGlkZWUxbdo0Tpw40e52VAsaAmgwS4RoALXI6/WydOlStm7dyr59+1i/fj379u0L\nW/0Oh4NHH32Uffv2sXPnTn7zm9+Etf7GVq9ezciRIzukboBly5YxY8YMPvnkE/73f/837G2Vlpby\n5JNPkp+fz969e/F6veTl5bW73ptvvplt27Y1eW/VqlVMmTKF/fv3M2XKlLD/4VSNNKRyBrNEiAbQ\nZtxuNx6Ph7ZSXP/xj38wbNgwhg4dSkxMDHPnzmXjxo1h60daWhoTJkwAICkpiZEjR1JaWhq2+huU\nlJSwefNmFi9eHPa6AU6ePMnf/vY3brnlFgBiYmLo3bt32NvxeDzU1tbi8Xioqalh4MCB7a7zW9/6\nFn37Nn0g5caNG1m4cCEACxcuZMOGDe1uR7UgzCNQY8wMY8ynxpgCY8zyM3x+jjHmbWPMHmPMP40x\n326rTg2gzXi9XgoLC9sMoqWlpWRmZgZeZ2RkdEiAAygqKmLPnj3k5OSEve477riDX/7yl9hsHfNT\nKCwsJCUlhe9+97uMHz+exYsXU11dHdY20tPT+clPfsI555xDWloavXr14oorrghrGw2OHDlCWloa\nAKmpqRw5cqRD2lF8dREpmKWtqoyxA78BrgRGAfOMMaOarXYf8IqIjAfmAk+1Va8G0DMoKSmhqKiI\nAwcOtDkS7WhVVVXMmTOHJ554gp49e4a17tdff53+/ftz4YUXhrXexjweD++//z5Llixhz549JCYm\nhv2w98SJE2zcuJHCwkLKysqorq5m3bp1YW3jTIwxGNOu52Wo1oR3BHoRUCAiB0WkHsgDrmm2jgAN\n/8l6AWVtVaoBtAUlJSUUFxdTUFCAz+f72ufp6ekUFxc3WT89PT2sfXC73cyZM4f58+cze/bssNYN\n8M4777Bp0yYGDx7M3Llzeeutt1iwYEFY28jIyCAjIyMwes7NzeX9998PaxtvvvkmQ4YMISUlBafT\nyezZs3n33XfD2kaDAQMGcPjwYQAOHz5M//79O6QdRagBNNkYk99oua1ZbelAcaPXJaffa2wFsMAY\nUwJsAX7YVhc1gLahpKSEAwcO4PE0nTd14sSJ7N+/n8LCQurr68nLy2PmzJlha1dEuOWWWxg5ciR3\n3nln2Opt7OGHHw6MtvPy8pg8eXLYR26pqalkZmby6aefArB9+3ZGjWp+5NQ+55xzDjt37qSmpgYR\nYfv27R12UWzmzJk8//zzADz//PNcc03zQYwKq+ADaLmIZDda1lpobR7wnIhkAN8GXjTGtBojNYAG\noaSkhMLCwiaH9A6Hgz59+jB9+nRGjhzJ9ddfz+jRo8PW5jvvvMOLL77I008/zbhx4xg3bhxbtmwJ\nW/2NzZgxo0PqbfDrX/+aiy66iPPPP58PPviAf/u3fwtr/Tk5OeTm5jJhwgR69uyJz+fjttuaD0BC\nN2/ePCZNmsSnn35KRkYGTz/9NMuXL+c//uM/yMrK4s0332T58q9di1DhEt6r8KVAZqPXGaffa+wW\n4BUAEdkBxAHJrXZRH6jclMvlavU+xauvnkV19ckQa3UC7ugsY5wgobZjsZyltqJ431kqA0lJfais\nPB5yuS6ofQ9UTjWSH+QZJfNo6w9UNsY4gM+AKfgD5y7gRhH5qNE6W4GXReQ5Y8xIYDuQLq0EySi/\nzz/6VFefBBPi3xEx4AixjMdAfIhlag0kh1im3ECahb+Lhw1khliu2EJbhy1uk5V9Z+U7CvW3AJw6\npReeghLGVE4R8RhjfgC8gX/M+oyIfGSMeQjIF5FNwF3A740xP8Y/WLy5teAJGkCVUtEqzKmcIrIF\n/8Whxu890Ojf+4BvhFKnBlClVHTqArnwUd49pVS3pQFUKaXaQZ9Ir5RSFugIVCmlLNIHKnc9R48e\nxeVyERsbq3nOSkWSjkC7nuTkZGw2G9XV1cTFxeFw6C5SKiI0gHY9xhhiYmJwOBy4XC7cbjdxcXE6\nGlWqszWkckYxDaAtsNlsJCQk4Ha7qa6uJjY2FqfTCTj9mUUhcfizVkItU2uhTLmFMoet/HFw+DOL\nOrwti9tkZd9Z+Y5C/i2APwVUtUlHoF2f0+lsMhoFd+elZfYJscwJi+mVwyykchYYGBtiuQ8ttFVg\ncZus7LvOSP8EC4G6mzL4H+cRxfRpTEEwxhAfH09MTEyku6LOEh39BKyzQheYE0lHoCHQC0oqXJpP\nVqfOoAscwusIVKkoodMon4FOa6yUCoZOo9xMFziE1wCqVJTQaZSbCfO0xh0hys8wKNW9detplDWV\nUykVLt1uGmW9iKSUao9uPY1yFziE1wCqVBTr1tMod4EAGuUD5Gjk7Ly0zBMWylhJryywmMr5YWe0\nZXGbrOy7Tkn/hJZSOefNm8df/vIXysvLycjI4MEHH2T58uVcf/31PP300wwaNIhXXnnFQntdmObC\nn23cnZeWOSTEMoUGxodYZo+Byy2kI75t4OoQy/0/C229bXGbrOy7zkj/hBYD9fr168/4/vbt20Nv\n42zQBc6BRnn3lFLdll6FV0opi3QEqpRSFmkAVUopizSAKqWUdaJX4ZVSKnRig/oof6CyBtBmXC4X\nXq8Xuz3K//QpdZYTAx57sLk+vg7tS0s0gDbj8Xioq6vD5/PhdDpxOp3YbJqwpVRnE2PwBv0Q8/oO\n7UtLNIA206NHDxISEhAR3G43tbW1AIFgCs7OyyoqtFBmT4hljMN/s3qojMN/Y3yHt2Vhm6zuu07J\nXgKdVC543ig/EtQA2oKG6Y1jYmLw+Xy43W5qamoANySHmH1SbnFiNCsZOFayg35sIZvmcQP/HmK5\ney209bjFbbKy76x8R6H+FsDCLKPdk2DwRnkupx6bBsFmsxEbG0tiYmKku6LOEjqpXNsEgwd7UEuk\naABVKgKCmVTu8ccfZ/To0YwZM4Z58+bhcrk6oWfRQzDUExvUEikaQJWKQqWlpTz55JPk5+ezd+9e\nvF4veXl5ke5Wp2o4hA9miRQNoEpFKY/HQ21tLR6Ph5qaGgYOHBjpLnW6cAZQY8wMY8ynxpgCY8zy\nFta53hizzxjzkTHmpbbq1ItISkWh9PR0fvKTn3DOOecQHx/PFVdcwRVXXBHpbnWqhnOg4WCMsQO/\nAaYBJcAuY8wmEdnXaJ0s4B7gGyJywhjT5uP/dQSqVBQ6ceIEGzdupLCwkLKyMqqrq1m3bl2ku9Wp\n/IfwjqCWIFwEFIjIQRGpB/KA5o/3vxX4jYicABCRL9uqVAOoUlHozTffZMiQIaSkpOB0Opk9ezbv\nvvtupLvVqfwXkWKCWoBkY0x+o+W2ZtWlA8WNXpecfq+x84DzjDHvGGN2GmPavFVCD+GVikLnnHMO\nO3fupKamhvj4eLZv3052dnaku9WpBEI5hC8XkfbuIAeQBVwGZAB/M8aMFZGK1goopaJMTk4Oubm5\nTJgwAYfDwfjx47nttuaDqrOdCfbwPBilQGaj1xmn32usBHhPRNxAoTHmM/wBdVeLPRQJKZPCQtpF\n1+JyudixY0eLn0+ePB1wh1irA/B0fBnjAAmxjM0BvlD7ZrGclTJWtqmz9relMgBO/KfhznrtSrka\nkZ0oT+ePDGrdb5rdu1sbgRpjHMBnwBT8gXMXcKOIfNRonRnAPBFZaIxJBvYA40TkWEv16gg0ZG5I\nC/HvyGEDw0IsU2BxAjYrqZIbLfxdvMbQq+5wSEVOxqaF3tY1FrfJyr6z8h2F+lsA/+9BBSVc93iK\niMcY8wPgDfxzfT4jIh8ZYx4C8kVk0+nPrjDG7AO8wE9bC56gAVQpFaXCnQsvIluALc3ee6DRvwW4\n8/QSFA2gSqmoJBjqonxaTg2gSqmo1BWexqQBVCkVlTSAKqVUO0TyUXXB0ACqlIpKEt77QDtEdPdO\nKdVt6SG8UkpZ5L8KHxPpbrRKA2gzx48fp6qqCmPMGRelVOfQQ/guqG/fvvTo0QMROeMCTguZJA5/\n1koorMxgaXP4s3BCYXf4s31C5XD4M4s6ui0r22R19s9QvyMcFrOKHMyYMSOoaT26Oz2E76JaHnG6\nrc3eODbEMh9anI3SwkyZoaZkgj8t87D0CqlMmjlpLf3TyuyfVvadle8o1N8CQLEJKnhWVFSwePFi\n9u7dizGGZ555hkmTJoXeXhel50CVUpYtW7aMGTNm8Oqrr1JfX396Wu3uQwOoUsqSkydP8re//Y3n\nnnsOgJiYGGJi2r6g8sADD9C3b1/uuOMOAO6991769+/PsmXLOrK7HaIrpHLqE+mVikKFhYWkpKTw\n3e9+l/Hjx7N48WKqq6vbLLdo0SJeeOEFAHw+H3l5eSxYsKCju9shdFZOpZQlHo+H999/nyVLlrBn\nzx4SExNZtWpVm+UGDx5Mv3792LNnD3/6058YP348/fr164Qed4xoD6B6CK9UFMrIyCAjI4OcnBwA\ncnNzgwqgAIsXL+a5557jiy++YNGiRR3ZzQ4Vzlk5O4qOQJWKQqmpqWRmZvLpp58CsH37dkaNGhVU\n2VmzZrFt2zZ27drF9OnTO7KbHSrMs3J2CB2BKhWlfv3rXzN//nzq6+sZOnQozz77bFDlYmJiuPzy\ny+nduzd2e3SP4NqiV+GVUpaMGzeO/Pz8kMv5fD527tzJH/7whw7oVedpmNY4mmkADZnTf2N8SBz+\nm65DYRz+m7tDYXP4byIPhZWMIvxJRWnmZIht2UNvy8o2Wdl3Vr4jHBZ+C4Bxhl4mSPv27eOqq65i\n1qxZZGVldVg7naErnAPVABoynVQOwKuTykXlpHKjRo3i4MGDHVZ/Z9JceKWUagc9B6qUUhZoKqdS\nSlmk50CVUsoi/1X46M6F1wCqlIpKegivlFLtoAFUKaUs0HOgSillUVe4D1QfJtIGr9dLfX09tbW1\nQT2PUalgzJgxI9JdiHoNqZzBLJES3eE9Aqqrq3G5XHi9XkQEu92O3W4nJiYGm82GTip3mk4qh+VJ\n5Ywz6AnlvF4v2dnZpKen8/rrr4feVhemh/BdkM1mw+FwEBsb2/Kkcskhpu+VW5h8rNjA+BDL7LE4\nmVqoqZLgD2pWJnuzkpZpZZus7Dsr31GovwXw/x6CtHr1akaOHEllZWXo7ZwF9BC+i4mPj8fhcOgc\n8CriSkpK2Lx5M4sXL450VyKiK0zpEd3hXalu7I477uCXv/wlp06dinRXIqIr3AeqI1ClotDrr79O\n//79ufDCCyPdlYjyYA9qiRQdgSoVhd555x02bdrEli1bcLlcVFZWsmDBAtatWxfprnUaH7aoT+XU\nEahSUejhhx+mpKSEoqIi8vLymDx5crcKng3CeQ7UGDPDGPOpMabAGLO8lfXmGGPEGJPdVp06AlVK\nRaVwngM1xtiB3wDTgBJglzFmk4jsa7ZeErAMeC+YenUEqlSUu+yyy7rdPaAAQljPgV4EFIjIQRGp\nB/KAa86w3v8BHgFcwVSqAVQpFaVCmtY42RiT32i5rVll6UBxo9clp9/7qjVjJgCZIrI52B7qIbxS\nKiqFeAhfLiJtnrNsiTHGBjwG3BxKOQ2gIXOGlEniZ2X2Roc/OyYUVmfyDDVVsqFcqLNlWk3LtDLD\nZqj7zup3FPJvAaDjZuU8mwiGuvDluZcCmY1eZ5x+r0ESMAb4y+kkmlRgkzFmpoi0OLe0BtCQuSE+\nxPS9WgN9QixzwsCQEMsUWkxhDHUGS/DnmltJsbQyW6aVbbKy76x8R6H+FsD/e1BtCvPTmHYBWcaY\nIfgD51zgxkBbIieB5IbXxpi/AD9pLXiCBlClVBQL11V4EfEYY34AvAHYgWdE5CNjzENAvohsslKv\nBlClVFQKdyqniGwBtjR774EW1r0smDo1gCqlopJg8PqiOxdeA6hSKiqJz1Dniu5UTg2gSqmoJGLw\nenQEqpRSoRM0gJ5NvF5vpLtVWySYAAARTUlEQVSgVLchYvC4ozuAaipnEDweD9XV1dTV1UW6K+os\n0dakcsXFxVx++eWMGjWK0aNHs3r16k7qWTQx+LyOoJZI0RFoKzweD3V1dRhjiIuLw26P7r+Gquto\na1I5h8PBo48+yoQJEzh16hQXXngh06ZNY9SoUZ3UwygggB7Cdy0iEgicNpuN+Pj407NxNnBayCRx\n+LNWQi1T2EnpnyHPYIm1FEurs2VaScu0su+sfEeWsoraTuVMS0sjLc0/62lSUhIjR46ktLS0ewVQ\nnwFXdIeo6O5dBJSXl+N2u88QOBu4wRFi+p7HQsqf1fRPKzNLDrOQjlhgYGyI5T600FaBxW3qjLTM\nWhP6bwH8v4cQFBUVsWfPHnJyckJvq6vzRLoDrdMA2kxKSgrx8fGR7oZSAFRVVTFnzhyeeOIJevbs\nGenudC7/A0GjmgZQpaKU2+1mzpw5zJ8/n9mzZ0e6O51PA6hSygoR4ZZbbmHkyJHceeedke5OZAjg\njnQnWqe3MSkVhd555x1efPFF3nrrLcaNG8e4cePYsmVL2wXPJgLUBblEiI5AlYpC3/zmNxGxcIHq\nbKKH8EopZZEGUKWUskgDqFJKWaQB9GzkDPlGaGsZKxYzY6xMjFZgJZvG4b8xvsPbsrhNnZJV5LDw\nWwCdVC4EGkDPNm4wIZ7cFwsZK1azl5JDLFNuIM3CxYrDFjOEQm3rsMVt6oysIo8J/bcA/t+DapsP\ncEW6E63TAKqUik56CK+UUhZpAFVKKYs0gCqlVDtoAFVKKQt0BKqUUhb5gNpId6J1+jCRIIkINTU1\nke6GOku0NScS+Kf9GD58OMOGDWPVqlWd0KsoI4A3yCVCNIAGwev1Ul1djdOpN0Cr8GhrTiSv18vS\npUvZunUr+/btY/369ezbt6+TehdFPEEuEaIBtA1ut5va2loSEhI0gKqwaWsE+o9//INhw4YxdOhQ\nYmJimDt3Lhs3buyk3kWJhnOgURxA9RxoK1wuFz6fj8TERIwxiAgJCT2pqQk1k8RK+qeVyeuc/iyc\nUMsctpiOGHKKpZW2LG6TlX1n5TuyklVknJSUlDBjxowWR6KlpaVkZmYGXmdkZPDee++F3lZXpheR\nuiafz0dtbS0Oh4OEhATAfw7UZrPx8sv/F7vdTmxsbIf2oaH9jh711tbWEhMT0+FTNtfW1hIbG9vC\nRH3hU1NTQ1xcXIe203A+3OpU1z/96U85dOgQLper1SDa7WkqZ9dz8uTJwH8Oh8O/e2w2G5WVlfh8\nvsBI1OPpuD+NDQ/S9fl81NfXd1g7DW34fL4ObQP821Rb2zmXVKurqzs8UIM/WFtpZ8WKFcyfP5/i\n4mKKi4vp1asXkyZNahJI09PTKS4uDrwuKSkhPT09LP3uUnQE2rXExsaSkJCAzWZDRLDb7YEglpCQ\n0OEjNa/XS11dHfHx8RjTsQ+d8Hg8eDwe4uLiOrQdgPr6eowxnXIeuaamhtjY2A7/rhqOVKyMRDds\n2MDPfvYzysvLcblcfPnll01GoxMnTmT//v0UFhaSnp5OXl4eL730UkdsRvTqAofwehGpmYbDv4bg\n6XK5cLlcxMfHd/h/SBHB5XIRFxfX4cET/AG0sy6MNYzcO0NsbCx1dR0/UY7NZiM+Ph6Xy4XXG/q9\nNI888ggpKSkkJCRQXFxMaWlp4OKSw+FgzZo1TJ8+nZEjR3L99dczevTocG9CdGuYVC6YJUJMiD/q\ns36SlpqaGnbs2IHdbsftdlNXV9dpAc3tdmOMCZw66Gh1dXUdfi63gc/nw+v1dlrArqurIyYmplO+\nt4ZTLVZH8suXL+fEiRNUVFSQlpZGenr62XJetF073yRnCzPzg1v5WbNbRLLb054VGkCbERH+/ve/\nWxpRKGXV3XffTUVFBbGxscTGxpKcnHw2BNH2BdB+2cK/BhlAX2w7gBpjZgCrATvwXyKyqtnndwKL\n8Z84OAosEpFDrdapAVQp1UHaF0D7ZgtTggygr7YeQI0xduAzYBpQAuwC5onIvkbrXA68JyI1xpgl\nwGUickNrzeo5UKVUdApvKudFQIGIHBSReiAPuKZJcyJvi0hDvvZOIKOtSvUqvFIqOoV2FT7ZGNN4\nuLpWRNY2ep0OFDd6XQLktFLfLcDWthrVAKqUik6hBdDycF1EMsYsALKBf2lrXT2EV6qDHT9+nGnT\nppGVlcW0adM4ceLEGdd7/vnnycrKIisri+effz7w/u7duxk7dizDhg3jRz/6UeB2sD/84Q+MHj0a\nm81Gfn7Tc4UPP/www4YNY/jw4bzxxhuB97vUE57CextTKZDZ6HXG6feaMMZMBe4FZopIm/fC6UUk\npTrY3XffTVxcHDt27GDPnj306NGDPXv20KdPn8A6x48fJzs7m7vuuovHHnuMzz//nNWrV3P77bdz\n0UUXcfvtt/Poo49y8OBBpk6dyoYNG/jkk0+orKxkxowZJCUlMXz4cF555RUOHz7MlClT6N+/P263\nm4KCArxeL1988QUXX3wxdXV19O7dmwMHDnDuueeyd+/ejtr09l1E6pktZAd5EentNi8iOfBfRJqC\nP3DuAm4UkY8arTMeeBWYISL7g2lWR6BKdbCNGzdSXl7OpEmTyMrK4tChQ/Tv35/LLrssMBp94403\nuPTSS3n00Ue58847cTgcLF26lOTkZIqLi3nqqae466676NGjB5s2baJHjx489dRTvPrqq/Ts2ZOU\nlBT++te/0r9/f6666ipuvfVWVq9eTVFRER6PBxFhzJgxDBgwAIfDwf333099fT2ffPIJAwcOJDu7\n02+hbFtDLnwwSxtExAP8AHgD+Bh4RUQ+MsY8ZIyZeXq1/wB6AH8wxnxgjNnUVr06AlWqg/Xu3ZsB\nAwaQlZXFW2+9FXg8YnZ2NhdffDGPPPIIv/rVr9ixYwe7d++mtLQUEeG8886jrKwMl8tFYmIi1dXV\neL1eUlJSsNvtDBgwgKKiIioqKrDb7VxyySXs3bsXp9PJiBEjOHz4MIWFhdjtdm688Ub69+/P2rVr\nqampob6+HhHh2muvZf/+/axfv55Ro0aFe9PbNwJNzBZGBTkCzY/MjfQ6AlUqDKZOncqYMWO+tjQ8\nw/OLL77gT3/6U+DRiDU1NXzwwQds2LAhUMfHH3+MzWajb9++GGPYv38/AwYMwG63U11dTUJCAiLC\nF198gdPppGfPnhw7dgybzYbX6+Wdd97h5MmTVFdXs2PHDhwOB3a7nZqaGt5++22GDBnCyZMnA89A\n8Hq9bNiwgZSUlOh81qg+kV6p7uHNN98kNTX1a+/fe++9JCYm4vF48Pl8HDt2jJSUFAAqKyspKioC\n/E9fKisrQ0RwOp0kJSXh8Xg4ePAgdXV1OBwOXC5XILgWFhZy8OBBwJ8C7HA46NmzJ263m6qqKlwu\nF/PmzQs8COfAgQP86Ec/AmD27NlkZGQQExOD1+vlr3/9K3/84x87YS9ZEOUPVNYAqlSYNATRoqIi\nCgoKAsuRI0eora0NXD13uVz07NkT8D+l6sc//jHTp0+nqqqK8vJybDYbVVVV2O32wEixtraWpKQk\njh8/TnJyMsYYDh06FMj1r6+vZ/r06QCBxxPeeuutJCcnY7PZsNls1NbW4na7effdd6mtrQ08ii8m\nJobdu3ezaNGizt5lresCT6TXc6BKdbBjx44xYsSIQHBMT0/n6NGjuN1uvF4vdrudH/7wh+Tl5fHF\nF19gjCEtLY2ysjLsdnuT5zLExMQEZkk4efJkk3ZsNhsxMTG4XP6rKgMHDqSsrIxZs2bx2muv4XQ6\n8Xg8GGPw+XwkJCRQU1PDoEGD+Pzzz0lJSWHr1q1MmDAhXJvevnOgsdlCepDnQAv1HKhSZ6V+/frx\n4osvBp5EdezYMeLi4hARkpOT6dOnD3l5eQwdOpR+/fohIpSVleFwOAKjyYaRptvtxul0UlVVRXx8\nfODJXcYYbDYbAwYMoEePHgCMGzeOjIwMNmzYEChvs9lISkoCCDw96vPPPychIYH58+ezZMmSztsx\nbekCI1ANoEp1gqlTp9KvXz9iY2NxuVxUVFTgdDq55JJLGDp0KC6Xi5KSEuLj4wOH6A23HzXcL2q3\n23E6ndTW1uLz+UhJSQmMOuPj4/F4PJSUlOB0OrHb7Rw8eJCjR48iIoEpaYwxXH311cTExARuoRIR\nvF4vOTk5VFRUcPjw4UjuqqY0gCqlHA4Ha9eupb6+HrvdzuDBg0lISAhcUc/OzubKK68kOTmZEydO\nEBMTQ0pKCldeeWVg1NhwKD9q1CiSk5MZPnw4F1xwAfX19dTV1TF+/HhGjx7NtddeS3x8PKWl/kSb\nWbNm0atXLy644ALGjh3LiRMnqK+vZ9CgQYwbN47+/fvTt29f7rnnHjIyMgLlIq4LPFBZA6hSneTq\nq6/mtddew+fzUVpayuDBg3G5XIGpsx0OBzt27GDs2LF4PB4qKyvJycnhsssuC9y25PF4OHDgANde\ney0FBQX07t2biRMnIiJ8/PHHzJ49mx07dhAbG0tVVRUDBgzg4MGDnDp1ivfff5/JkyczY8YMjDEc\nO3aMwsJCFi5cSHp6Oj6fD7c7gtGoOb2NSSnVWOMg2hDwGoJoZWVlIOWzIYi+/PLLXHrppYEgevvt\nt+PxePjjH//I/fffHwiiq1atwuPx8PTTTweu6H/44Yece+65eL3ewIWpFStWMGfOHLKyspgzZw4+\nn4/169ezdOlSfD4fR44ciZ7J67rAOdDA+ZEgF6VUGGzatEkcDodkZmbKihUr5Pzzz5fvfe97snHj\nRhERqa2tlXHjxknPnj1l4sSJ8sQTT8h1110nIiI//OEPJSYmRoYNGybPPPOMDBkyRDwej2zevFmy\nsrIkMzNT+vfv36S9jIwMmTRpklxwwQWSmpoqPXv2lGHDhklmZqYkJSXJRRddJP/5n/8pEydODOdm\nhhpfmixwoeCQ4BbIb297lvoYYgGlVJg0BLyhQ4fKypUrRUTk/vvvbxJEc3Nz5dxzz5WJEyfKgQMH\nAmVXrlwpQ4cOlfPOO0+2bNkSeH/u3LmSmpoqDodD0tPT5b/+679ERKS8vFwmT54sw4YNk8mTJ0tO\nTo5s2rRJfD6f3H777TJ06FAZM2aM7Nq1K5yb2P4AaiS4JUIBVO8DVaqbuf3220lJSeHBBx/s6Kba\ndx+oyRYI8j5QInMfqD5QWalu5LnnnuPQoUOsWbMm0l05K2gAVaqb2L17N7/61a/4n//5n0Aap2of\n3YtKdRNr1qzh+PHjXH755YwbN47FixdHuktdnp4DVUp1lHaeA50g8E6QayfoOVCllPpKQypS9NIA\nqpSKUqFNyxkJGkCVUlFKR6BKKWWRBlCllLJIgNpId6JVGkCVUlFKz4EqpZRFegivlFIW6QhUKaUs\n0hGoUkpZpCNQpZSyyIdehVdKKUv0EF4ppdpBD+GVUsoCHYEqpZRFGkCVUsoivQqvlFIW6VV4pZSy\nSA/hlVLKoug/hNdJ5ZRSUaphBBrM0jZjzAxjzKfGmAJjzPIzfB5rjHn59OfvGWMGt1WnBlClVJRq\nGIEGs7TOGGMHfgNcCYwC5hljRjVb7RbghIgMAx4HHmmrXg2gSqko1XARKZilTRcBBSJyUETqgTzg\nmmbrXAM8f/rfrwJTjDGtziwa6jnQdk1TqpRSwTv8BqxIDnLlOGNMfqPXa0VkbaPX6UBxo9clQE6z\nOgLriIjHGHMS6AeUt9SoXkRSSkUlEZkR6T60RQ/hlVLdQSmQ2eh1xun3zriOMcYB9AKOtVapBlCl\nVHewC8gyxgwxxsQAc4FNzdbZBCw8/e9c4C0RkdYq1UN4pdRZ7/Q5zR8AbwB24BkR+cgY8xCQLyKb\ngKeBF40xBcBx/EG2VaaNAKuUUqoFegivlFIWaQBVSimLNIAqpZRFGkCVUsoiDaBKKWWRBlCllLJI\nA6hSSln0/wG61jFFXnh5ZAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1a0ff57710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solving time step:  3\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVAAAADxCAYAAACd3+8mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl0VFW6+P3vrqrMhDGBQBIZJDKq\ngMGItt3KINityBAVBEUR7Ua08aqteHEAX1qxr6gobXuxHeGVaNtXYMmgLdr2VUGJ4m0RB4YEMgAS\nIIQMldTw/P4IKZNIkqqTSqoCz2ets1aq6py996mqPLXP8OxtRASllFKBs4W6AUop1VZpAFVKKYs0\ngCqllEUaQJVSyiINoEopZZEGUKWUskgDqFLqlGeMeckY86MxZnsDrxtjzDPGmF3GmH8bY4b5U64G\nUKXU6eAVYFwjr18OpJ1YbgX+4k+hGkCVUqc8EfkXcKSRVa4CXpNqW4COxpjuTZXrCLQdAa6vlDp9\nmeZs3NcYKfdz3f3wDeCs9dRyEVkeQHXJQF6tx/knntvf2EaBBlCllGoV5cBv/Vx3AThFJL0Fm3NS\nGkCVUmHJ0KoBqgBIrfU45cRzjdJzoGGmXbt27Nmz56SvvfLKK/ziF79o5RYpFRo2IMbPJQjWAjec\nuBp/AXBMRBo9fK9p42nn9ddfJz09nXbt2tG9e3cuv/xyPv74Y8vlGWPYtWtXneeOHz/OXXfdRa9e\nvYiLi+OMM84gMzOTzz77rNGySktL6dOnj6V2VFVVsWDBAtLS0oiLi6NXr17MnDmT3NxcS+W1lAUL\nFjB9+vRQN6NF5ebmYozB7XaHuiltlgEi/FyaLMuYVcBmoJ8xJt8Yc7Mx5nfGmN+dWGU9sAfYBbwA\n3OZPG0+7Q/gnn3ySxYsX8/zzzzN27FgiIyPZuHEja9asCbh353a7cTh+/hZWVlYycuRIOnbsyDvv\nvMOAAQNwOp1s2LCBDRs2kJGR4XdZgcjMzCQ/P5/XX3+doUOHUlZWxsqVK9m0aRM333xzs8pWwReM\nz/xUFsxDeBGZ2sTrAsyxUnAgS5tWXFwscXFx8uabbza4zmeffSYXXHCBdOjQQZKSkmTOnDlSWVnp\nex2QZcuWSd++faVXr15y8cUXCyCxsbESFxcnWVlZ8sILL0hSUpKUlpY22p76ZdU8t3PnThERKSoq\nkiuvvFLi4+Nl+PDh8sADD8hFF1100rL+8Y9/SHR0tOzbt6/B+goKCuTKK6+UTp06yZlnninLly/3\nvfbwww9LZmamTJs2Tdq1ayeDBw+W77//Xh599FFJTEyUlJQUeffdd33r/+pXv5J58+bJ8OHDJT4+\nXsaPHy+HDx8WEZEPP/xQkpOT69Tds2dP+cc//iEbNmyQiIgIcTgcEhcXJ+ecc46IVH82M2fOlKSk\nJOnRo4fMnz9f3G73SfejvLxcbrjhBunYsaP0799fHn/88Tr1FRQUyKRJkyQhIUF69eolS5cu9b3m\ndDpl7ty50r17d+nevbvMnTtXnE5nnXY//vjjkpiYKElJSfL222/LunXrJC0tTTp16iR//OMffWV5\nPB557LHHpE+fPtK5c2e5+uqrfe9BamqqABIXFydxcXHy6aefyssvvywXXnih3HnnndK5c2e5//77\npVOnTvLvf//bV+bBgwclJiZGfvzxxwY/xzYk0PhSZ+kJ8oKfC5Dd3PqsLKdVAN2wYYPY7XZxuVwN\nrpOdnS2bN28Wl8slOTk50r9/f3nqqad8rwMyevRoOXz4sJSXl/ueqwl6IiLXXnutzJgxo8n2NFXW\ntddeK1dffbWUlpbK119/LT169GgwgN53333yy1/+stH6Lr74Ypk9e7ZUVFTItm3bJCEhQTZt2iQi\n1QE0KipKNm7cKC6XS66//nrp1auXLFq0SKqqqmT58uW+IC9SHUB79OghX3/9tZSWlsqkSZNk2rRp\nItJ4AK2pq2bdGhMmTJBbb71VSktL5eDBgzJ8+HB5/vnnG93XI0eOSF5enpx99tm++jwejwwbNkwW\nLlwolZWVsnv3bundu7ds3LhRREQefPBBycjIkIMHD8qPP/4oI0aMkAceeMDXbrvdLgsXLvTtc0JC\ngkydOlVKSkpk+/btEh0dLXv27BERkaeffloyMjIkLy9PnE6n3HrrrTJlyhQREcnJyRGgznft5Zdf\nFrvdLs8884y4XC4pLy+X2bNny7333utb5+mnn5Yrrrii0c+xDWlWcOoF8rKfiwbQVrBy5Urp1q1b\nQNs89dRTMmHCBN9jwBd0aj9XO4COGjVK7rvvPt/jbdu2SYcOHSQ+Pl7OOussv8pyu93icDjk22+/\n9b12//33NxhAZ82aJddee22D+7Fv3z6x2WxSUlLie27evHm+QP/www/L6NGjfa+tXbtW4uLifL3A\nkpISAeTo0aMiUh1Aa+/jN998IxEREeJ2uwMOoAcOHJDIyEjfj4iIyOuvvy6XXHLJSfeldkAUEXnh\nhRd89W3ZskVSU1PrrP/oo4/KjTfeKCIiffr0kXXr1vle27hxo/Ts2VNEqgNodHT0z/Z5y5YtvvWH\nDRsmb7/9toiI9O/fX95//33fa4WFheJwOHw/vicLoPXbVtNer9crIiLnnXeevPHGGyfd7zaoWcGp\nN8hKP5dQBdDT6gRMly5dKCoqavTc0w8//MBdd91FdnY25eXluN1uzjvvvDrrpKamnnTb2vXs3//T\nBbwhQ4ZQXFzM+++/z6xZs/wq69ChQ7jd7jqv9+zZs9E6f/jhhwZfLywspHPnzsTHx9cpLzs72/e4\nW7duvr9jYmJISEjAbrf7HkP1Ra6OHTv+rO09e/bE5XJRVFTUYBsasnfvXlwuF927/5T44fV6G3xv\nCgsL67xW+++9e/dSWFjoayOAx+Ph4osv9m1b+33s2bMnhYWFvsddunT52T7Xf19KS0t9dU2cOBGb\n7adrsXa7nYMHDza4r/X3KSMjg9jYWP75z3/SvXt3du3axfjx4xvc/nRScxU+nJ1WV+FHjBhBVFQU\nq1evbnCd2bNn079/f3bu3ElJSQmPPvpodVe9FmMaT7AYNWoU7733HmVlZU22qaGyEhMTcTgc5OX9\nlByxb9++BssZPXo0n3/+Ofn5+Sd9vUePHhw5coTjx4/XKS85ObnJNjakftsiIiJISEggLi6O8vKf\nckg8Hg+HDh3yPa6/z6mpqURFRVFUVERxcTHFxcWUlJTwzTffnLTe7t2719nP2u1ITU2ld+/evnKK\ni4s5fvw469evB6rfh71799Zpd48ePSztf2pqKhs2bKhTl9PpJDk5ucHP9WTPz5gxg5UrV7JixQoy\nMzOJjo621J5TTTCvwreU0yqAdujQgUceeYQ5c+awevVqysvLcblcbNiwgXvvvReovv2offv2tGvX\nju+++46//KXpMQW6detW597NG264ge7duzNx4kS2b9+Ox+PB6XTW6e01xW63M2nSJBYsWEB5eTk7\nduzg1VdfbXD90aNHM2bMGCZOnMgXX3yB2+3m+PHjPP/887z00kukpqZy4YUXcv/99+N0Ovn3v//N\niy++2KzbiVauXMmOHTsoLy/noYceIjMzE7vdzllnnYXT6WTdunW4XC4WLVpEZWWlb7tu3bqRm5uL\n1+sFqgPiZZddxt13301JSQler5fdu3fz0UcfnbTea665hscee4yjR49SUFDAsmXLfK+df/75xMfH\n8/jjj1NRUYHH42H79u1s3boVgKlTp7Jo0SIOHTpEUVERjzzyiOX34He/+x3z58/3BeRDhw6xZs0a\noPoH0GazNXhPb23Tp0/n7bffZuXKldxwww2W2nKqcvi5hMppFUAB7r77bp588kkWLVpEYmIiqamp\nLFu2jAkTJgDwxBNP8PrrrxMfH88tt9zCtdde22SZCxYsYMaMGXTs2JE333yT6OhoPvzwQwYOHMhv\nfvMb2rdvT79+/di6dStvvvmm321dtmwZpaWlJCUlceONN3LTTTc1uv5bb73Fr3/9a6699lo6dOjA\n4MGDyc7OZvTo0QCsWrWK3NxcevTowcSJE1m4cKHvNSuuv/56brzxRpKSknA6nTzzzDNA9Q/Vc889\nx6xZs0hOTiYuLo6UlBTfdldffTVQfbg8bFj1qGGvvfYaVVVVDBw4kE6dOpGZmVnnNEhtDz30ECkp\nKfTu3ZvRo0eTmZlJVFQUUP3D88477/DVV1/Ru3dvEhISmDVrFseOHQPggQceID09nXPOOYezzz6b\nYcOG8cADD1ja/7lz5zJ+/Hguu+wy4uPjueCCC3z3+cbGxjJ//nwuuugiOnbsyJYtWxosJzU1lWHD\nhmGM8Z1qUG2jB2rqH542QQcTUQBccsklTJ8+/WfndEPhL3/5C1lZWQ32WNuCmTNn0qNHDxYtWhTq\npgRTswYT6WeM/Lef614KX4jmwivVtP3797Nnzx5GjBjBzp07WbJkCbfffnuom2VZbm4u//M//8O2\nbdtC3ZSwoheRlGoBVVVV/Pa3vyU+Pp6RI0dy1VVXcdttfmXehZ0HH3yQwYMH84c//IHevXuHujlh\nRQ/hlVKns2Ydwg8yRlb5ue65egivlFI/qemBhjMNoEqpsNTK44FaEu7tU0qdprQHqpRSFhnC/yq8\nBlClVFgyQIS/ESpE41ZrAFVKhSVjwO/xpjWAKqXUT4yBCHuoW9E4vZFeBcXWrVs555xzcDqdlJWV\nMWjQILZv3x7qZqk2rKYH6s8SsjbqjfQqWB544AGcTicVFRWkpKRw//33h7pJKrSadSN9eoSR7M5+\nVvRjaG6k1wCqgqaqqorhw4cTHR3Np59+6huYWJ22mhdAI41kJ/pZUaFmIqk27vDhw5SWluJyuXA6\nncTFxYW6SaotawN30msPVAXN+PHjmTJlCjk5Oezfv7/OQMfqtNS8HmiUkeyUptcDMHu0B6rasNde\ne42IiAiuu+46PB4PF154IR988AEjR44MddNUW2WAMD8LpD1QpVRLaV4PNMZItp8j/JlvtQeqlFI/\nMUBUqBvROA2gSqnw1AYuIoV585RSpy0NoEop1QxhfhFJA6hSKjxpD1QppSzSAKqUUhbpVXillLKo\nDfRAdTi7Zti4cSP9+vWjb9++LF68OKhl5+XlcemllzJw4EAGDRrE0qVLg1p+bR6Ph6FDh3LFFVe0\nSPnFxcVkZmbSv39/BgwYwObNm4Nex1NPPcWgQYMYPHgwU6dOxel0NrvMmTNn0rVrVwYPHux77siR\nI4wZM4a0tDTGjBnD0aNHm12PakBNAPVnCRENoBZ5PB7mzJnDhg0b2LFjB6tWrWLHjh1BK9/hcLBk\nyRJ27NjBli1b+POf/xzU8mtbunQpAwYMaJGyAebOncu4ceP47rvv+L//+7+g11VQUMAzzzxDdnY2\n27dvx+PxkJWV1exyb7zxRjZu3FjnucWLFzNq1Ch27tzJqFGjgv7DqWqpSeX0ZwkRDaD1uFwu3G43\nTaW4fv755/Tt25c+ffoQGRnJlClTWLNmTdDa0b17d4YNGwZAfHw8AwYMoKCgIGjl18jPz2fdunXM\nmjUr6GUDHDt2jH/961/cfPPNAERGRtKxY8eg1+N2u6moqMDtdlNeXk6PHj2aXeYvf/lLOneuOyDl\nmjVrmDFjBgAzZsxg9erVza5HNSDIPVBjzDhjzPfGmF3GmHknef0MY8yHxphtxph/G2N+3VSZGkDr\n8Xg85OTkNBlECwoKSE1N9T1OSUlpkQAHkJuby7Zt28jIyAh62XfeeSd/+tOfsNla5quQk5NDYmIi\nN910E0OHDmXWrFmUlZUFtY7k5GTuuecezjjjDLp3706HDh247LLLglpHjYMHD9K9e3cAkpKSOHjw\nYIvUo/jpIpI/S1NFGWMH/gxcDgwEphpjBtZb7QHgTREZCkwBnmuqXA2gJ5Gfn09ubi67d+9usifa\n0kpLS5k8eTJPP/007du3D2rZ77zzDl27duW8884Larm1ud1uvvzyS2bPns22bduIi4sL+mHv0aNH\nWbNmDTk5ORQWFlJWVsbKlSuDWsfJGGMwplnjZajGBLcHej6wS0T2iEgVkAVcVW8dAWr+yToAhU0V\nqgG0Afn5+eTl5bFr1y68Xu/PXk9OTiYvL6/O+snJyUFtg8vlYvLkyUybNo1JkyYFtWyATz75hLVr\n19KrVy+mTJnCBx98wPTp04NaR0pKCikpKb7ec2ZmJl9++WVQ63j//ffp3bs3iYmJREREMGnSJD79\n9NOg1lGjW7du7N+/H4D9+/fTtWvXFqlHEWgATTDGZNdabq1XWjKQV+tx/onnalsATDfG5APrgTua\naqIG0Cbk5+eze/du3O6686YOHz6cnTt3kpOTQ1VVFVlZWYwfPz5o9YoIN998MwMGDOCuu+4KWrm1\nPfbYY77edlZWFiNHjgx6zy0pKYnU1FS+//57ADZt2sTAgfWPnJrnjDPOYMuWLZSXlyMibNq0qcUu\nio0fP55XX30VgFdffZWrrqrfiVFB5X8ALRKR9FrLcgu1TQVeEZEU4NfACmNMozFSA6gf8vPzycnJ\nqXNI73A46NSpE2PHjmXAgAFcc801DBo0KGh1fvLJJ6xYsYIXX3yRIUOGMGTIENavXx+08msbN25c\ni5Rb49lnn+X888/nnHPO4auvvuI///M/g1p+RkYGmZmZDBs2jPbt2+P1ern11vodkMBNnTqVESNG\n8P3335OSksKLL77IvHnz+K//+i/S0tJ4//33mTfvZ9ciVLAE9yp8AZBa63HKiedquxl4E0BENgPR\nQEKjTdQBletyOp2N3qd45VUTKTt+LLBCbRHgdZ0627RmXafaNkB8+06UHDsS8HZtUPMGVE4yku3n\nGSWzpPEBlY0xDuAHYBTVgXMrcJ2IfFNrnQ3AGyLyijFmALAJSJZGgmSY3+cffsqOH4OpAf6OrDJw\nfYDbrDAwK8Bt/mrgjgC3edbAPRZ+F58wcH+A2z1moa4nLO6TlffOymcU6HcBOL5KLzz5JYipnCLi\nNsbcDrxLdZ/1JRH5xhjzCJAtImuBu4EXjDH/QXVn8cbGgidoAFVKhasgp3KKyHqqLw7Vfu6hWn/v\nAC4KpEwNoEqp8NQGcuHDvHlKqdOWBlCllGoGHZFeKaUs0B6oUkpZpAMqtz2HDh3C6XQSFRWlec5K\nhZL2QNuehIQEbDYbZWVlREdH43DoW6RUSGgAbXuMMURGRuJwOHA6nbhcLqKjo7U3qlRrq0nlDGOa\nyllP/VROl8tFZWUlUVFRREREMHL02MDT94wDxN30eqHYxuYAb4DbWN3Oyjbh/N5Z2QbAFoF4qgLf\nru1pXipnHyPZj/hZ0fWNp3K2FO2BNiEiIqJObxSvy1rK300BbvOyxRRGK+mVCy38Lj5s4I8Bbjff\nQl0PW9wnK++dlc8o0O8CVH8fVNMM1cN5hDEdjckPxhhiYmKIjIwMdVPUKaKlR8A6JbSBOZG0BxoA\nvaCkgqX+ZHXqJNrARSTtgSoVJnQa5ZPQaY2VUv7QaZTraQOH8BpAlQoTOo1yPUGe1rglhPkZBqVO\nb6f1NMqayqmUCpbTbhplvYiklGqO03oa5TZwCK8BVKkwdlpPo9wGAqimctbT1KycYZ3K2Vrpla1Z\n16mW/gkNpnJOnTqVf/7znxQVFdGtWzcWLlzIhAkTuOaaa9i3bx89e/bkzTff/NmFpjDWvFTOQUay\n3/CzorM1lbNt8Lpab7ZMKymMSwLc5m4DL1r4XbzZQFaA202xUNfNFveptdI/A/0uQPX34SRWrVp1\n0uc3bdoUeB2ngjZwDjTMm6eUOm3pVXillLJIe6BKKWWRBlCllLJIA6hSSlknYT4ivQZQpVRYEhtU\nhfmAyhpA63E6nXg8Huz2MP/pU+oUJwbcdn9zfbwt2paGaACtx+12U1lZidfrJSIigoiICGw2TdhS\nqrWJMXj8HsQ8NHNMaSZSPTWZSCKCy+WqngcJfMF01JhxrZOJFM6ZPgB2B3jCtH3hnL0EOqmcn4am\n2+SDbP+O4TubCs1ECic10xtHRkbi9XpxuVyUl5dXB89wziqykh30hYXfxfMMFDgD2yY5OvC6zrO4\nT+GavQTV3wfVJMHgCfN5jfXY1A82m42oqCji4uJC3RR1itBJ5ZomGNzY/VpCRQOoUiHgz6RyTz31\nFIMGDWLw4MFMnToVpzPAHn8bJxiqiPJrCRUNoEqFoYKCAp555hmys7PZvn07Ho+HrKysUDerVdUc\nwvuzhIoGUKXClNvtpqKiArfbTXl5OT169Ah1k1pdMAOoMWacMeZ7Y8wuY8y8Bta5xhizwxjzjTHm\n9abK1ItISoWh5ORk7rnnHs444wxiYmK47LLLuOyyy0LdrFZVcw40GIwxduDPwBggH9hqjFkrIjtq\nrZMG3A9cJCJHjTFNDv+vPVClwtDRo0dZs2YNOTk5FBYWUlZWxsqVK0PdrFZVfQjv8Gvxw/nALhHZ\nIyJVQBZQf3j/W4A/i8hRABH5salCNYAqFYbef/99evfuTWJiIhEREUyaNIlPP/001M1qVdUXkSL9\nWoAEY0x2reXWesUlA3m1HuefeK62s4CzjDGfGGO2GGOavFVCD+GVCkNnnHEGW7Zsoby8nJiYGDZt\n2kR6eqvfJx5SAoEcwhcF4UZ6B5AGXAKkAP8yxpwtIsWNbaCUCjMZGRlkZmYybNgwHA4HQ4cO5dZb\n63eqTnXG38NzfxQAqbUep5x4rrZ84DMRcQE5xpgfqA6oWxtsoaZy1tUik8q1VjqilfRKK9sAOBzg\nDtP2hXP6J2gqp5/6p8fJi9kD/Fr3F+aLRlM5jTEO4AdgFNWBcytwnYh8U2udccBUEZlhjEkAtgFD\nRORwQ+VqDzRQXhfcE+DvyBMGFga4zcMWJ2CzkCp5pmwPbBtgtxnM9bI8oG1WmFsDrmu3GWwt/dPK\ne2flMwr0uwDV3wfll2Dd4ykibmPM7cC7gB14SUS+McY8AmSLyNoTr11mjNkBeIA/NBY8QQOoUipM\nBTsXXkTWA+vrPfdQrb8FuOvE4hcNoEqpsCQYKsN8Wk4NoEqpsNQWRmPSAKqUCksaQJVSqhlCOVSd\nPzSAKqXCkgT3PtAWEd6tU0qdtvQQXimlLKq+Ch8Z6mY0SgNoPUeOHKG0tBRjzEkXpVTraAuH8JrK\nWU/tWTlPtoz79ZWnViqnww5uT2DbAMZhQ9wBzsVtpa5wTuW0OqOpzcHYMaP8mtajjWtWj6NXeheZ\nn/0bv9a91azQWTnDSYM9Tq/L2uyNfwxwm/kWZ6O0MFNmoCmZUJ2WKWsD28aM91hK/7Q0+6eV987K\nZxTodwHgMeNX8CwuLmbWrFls374dYwwvvfQSI0aMCLy+NkrPgSqlLJs7dy7jxo3jrbfeoqqqqnpa\n7dOIBlCllCXHjh3jX//6F6+88goAkZGRREY2fUHloYceonPnztx5550AzJ8/n65duzJ37tyWbG6L\naAupnDoivVJhKCcnh8TERG666SaGDh3KrFmzKCsra3K7mTNn8tprrwHg9XrJyspi+vTpLd3cFqGz\nciqlLHG73Xz55ZfMnj2bbdu2ERcXx+LFi5vcrlevXnTp0oVt27bx3nvvMXToULp06dIKLW4Z4R5A\n9RBeqTCUkpJCSkoKGRkZAGRmZvoVQAFmzZrFK6+8woEDB5g5c2ZLNrNFBXNWzpaiPVClwlBSUhKp\nqal8//33AGzatImBAwf6te3EiRPZuHEjW7duZezYsS3ZzBYV5Fk5W4T2QJUKU88++yzTpk2jqqqK\nPn368PLLL/u1XWRkJJdeeikdO3bEbg/vHlxT9Cq8UsqSIUOGkJ2dHfB2Xq+XLVu28Le//a0FWtV6\naqY1DmcaQANli6i+MT6gbRzVN10Hwu6ovrk7EA5H9U3kATAOW/XN6gFy2MGMD2wbS3VZ2CdL752V\nz8jmCPy7ANXfoRayY8cOrrjiCiZOnEhaWlqL1dMa2sI5UA2ggTrFJpUTnVSu2ikyqdzAgQPZs2dP\ni5XfmtpCLnx4t04pdVrTc6BKKWWBpnIqpZRFeg5UKaUsqr4KH9658BpAlVJhSQ/hlVKqGTSAKqWU\nBXoOVCmlLGoL94HqYCJN8Hg8VFVVUVFR4dd4jEr5Y9y4caFuQtirSeX0ZwkVnVSunsOHD/P555/j\n8XgQEex2u2+x2WyMGjNOJ5VDJ5WzvA2ALQLxVPm1qsfjIT09neTkZN55553A6wqtZqVcdUjvKxdl\n/8mvdTeYyTqpXDiw2Ww4HA6ioqIanlTujgB/R561MPnYYwaWBLjN3RYnogs0VZLqFFBLk71ZScu0\nsk9W3jsrn1Gg3wWo/j74aenSpQwYMICSkpLA6zkF6CF8GxMTE4PD4dA54FXI5efns27dOmbNmhXq\npoREW5jSI7zDu1KnsTvvvJM//elPHD9+PNRNCYm2cB+o9kCVCkPvvPMOXbt25bzzzgt1U0LKjd2v\nJVS0B6pUGPrkk09Yu3Yt69evx+l0UlJSwvTp01m5cmWom9ZqvNjCPpVTe6BKhaHHHnuM/Px8cnNz\nycrKYuTIkadV8KwRzHOgxphxxpjvjTG7jDHzGllvsjFGjDFNXtXXHqhSKiwF8xyoMcYO/BkYA+QD\nW40xa0VkR7314oG5wGf+lKs9UKXC3CWXXNIW7wFtNiGo50DPB3aJyB4RqQKygKtOst7/BzwO+HWP\nngZQpVSYCmha4wRjTHatpf7kW8lAXq3H+See+6k2Y4YBqSKyzt8W6iG8UiosBXgIX9ScTCRjjA14\nErgxkO00gAbKFhFQJkn1NhZmb7Q5qrNjAmFlNkq7ozrbJ1BWZ8sMtC6rM2wG+t5Z+YyMI/DvArTo\nrJynEsFQGbw89wIgtdbjlBPP1YgHBgP/PJFEkwSsNcaMF5EG55bWABoorwtmBZi+91cLKX+tmf4Z\n6AyWUD2LpZUUSyuzZYZrWuazJvDvAlR/H1STgjwa01YgzRjTm+rAOQW4zleXyDEgoeaxMeafwD2N\nBU/QAKqUCmPBugovIm5jzO3Au4AdeElEvjHGPAJki8haK+VqAFVKhaVgp3KKyHpgfb3nHmpg3Uv8\nKVMDqFIqLAkGjze8c+E1gCqlwpJ4DZXO8E7l1ACqlApLIgaPW3ugSikVOEED6KnE4wl86gullDUi\nBrcrvAOopnL6we12U1ZWRmXNuTTTAAARQklEQVRlZaibok4RTU0ql5eXx6WXXsrAgQMZNGgQS5cu\nbaWWhROD1+PwawkV7YE2wu12U1lZiTGG6Oho7Pbw/jVUbcfGjRsbfd3hcLBkyRKGDRvG8ePHOe+8\n8xgzZgwDBw5spRaGAQHC/BBeZ+Wsp6Kigo8//pjKykpsNhtRUVHYbD911EeOHhv4rJzGARKms0Ra\nnVnyVJst08pnZGUbCGhWzhpXXXUVt99+O2PGjAm8vtBpVsqVGZwu/E+jiUA/6Wd0Vs5wUFRUhMvl\nIiYmpk7g9PG64PoAf0dWGLgpwG1ebsX0z4UWfhcfNvDHALebb6Guh1sxLdPKZxTodwGqvw8ByM3N\nZdu2bWRkZAReV1tn4fepNWkArScxMZGYmJhQN0MpAEpLS5k8eTJPP/007du3D3VzWlf1gKBhTQOo\nUmHK5XIxefJkpk2bxqRJk0LdnNanAVQpZYWIcPPNNzNgwADuuuuuUDcnNAQI8HJDa9PbmJQKQ598\n8gkrVqzggw8+YMiQIQwZMoT169c3veGpRIBKP5cQ0R6oUmHoF7/4BQHeIXPq0UN4pZSySAOoUkpZ\npAFUKaUsagMBVDOR6nE6nWzevLnB11stEymcs5das65wzipqxUykNqp5mUh904U/+ZmJNFkzkdoG\nrwumBvg7sspCxsoKCxOWWZ287h4Lv4tPWMwQCrSuJ1ppsre/WvyMAv0uQPX3QTXNCzhD3YjGaQBV\nSoWnNnAIrwFUKRWeNIAqpZRFGkCVUqoZNIAqpZQF2gNVSimLvEBFqBvROB1MxE8iQnl5eaiboU4R\nTc2JBNXTfvTr14++ffuyePHiVmhVmBHA4+cSIhpA/eDxeCgrKyMiIiLUTVGniKbmRPJ4PMyZM4cN\nGzawY8cOVq1axY4dO1qpdWHE7ecSIhpAm+ByuaioqCA2NlYDqAqapnqgn3/+OX379qVPnz5ERkYy\nZcoU1qxZ00qtCxM150DDOIDqOdBGOJ1OvF4vcXFxGGMQEWLj2lMeaCaJLSLgeXCwRVRnxwS6zbMW\ntnnCQmaMLaI6s6il67K6T1beOyufkZWsIlsE+fn5jBs3rsGeaEFBAampqb7HKSkpfPbZZ4HX1Zbp\nRaS2yev1UlFRgcPhIDY2Fqg+B2qz2Xgj6//HbrcTFRXVom2oqb+le70VFRVERka2+JTNFRUVP5vh\ntCWUl5cTHR3dovXUnA+3OtX1H/7wB/bu3YvT6Ww0iJ72NJWz7Tl27Jjvn8PhqH57bDYbJSUleL1e\nX0/U7W65n8aaAV68Xi9VVS076ITX68Xr9bZoHVC9TxUVrXNJtaysrMUDNVQHayv1LFiwgGnTppGX\nl0deXh4dOnRgxIgRdQJpcnIyeXl5vsf5+fkkJycHpd1tivZA25aoqChiY2Ox2WyICHa73RfEYmNj\nW7yn5vF4qKysJCYmBmNadtAJt9uN2+0mOjq6ResBqKqqwhjTKueRy8vLiYqKavHPquZIxUpPdPXq\n1dx3330UFRXhdDr58ccf6/RGhw8fzs6dO8nJySE5OZmsrCxef/31ltiN8NUGDuH1IlI9NYd/NcHT\n6XTidDqJiYlp8X9IEcHpdBIdHd3iwROqA2hrXRir6bm3hqioKCorW36iHJvNRkxMDE6nE48n8Htp\nHn/8cRITE4mNjSUvL4+CggLfxSWHw8GyZcsYO3YsAwYM4JprrmHQoEHB3oXwVjOpnD9LiOh4oPWU\nl5ezefNm7HY7LpeLysrKVgtoLpcLY4zv1EFLq6ysbPFzuTW8Xi8ej6fVAnZlZSWRkZGt8rnVnGqx\n2pOfN28eR48epbi4mO7du5OcnHyqnBdt3nigCenCeD/HA305NOOBagCtR0T4+OOPLfUolLLq3nvv\npbi4mKioKKKiokhISDgVgmjzAmiXdOE3fgbQFU0HUGPMOGApYAf+KiKL671+FzCL6hMHh4CZIrK3\n0TI1gCqlWkjzAmjndGGUnwH0rcYDqDHGDvwAjAHyga3AVBHZUWudS4HPRKTcGDMbuERErm2sWj0H\nqpQKT8FN5Twf2CUie0SkCsgCrqpTnciHIlKTr70FSGmqUL0Kr5QKT4FdhU8wxtTuri4XkeW1HicD\nebUe5wMZjZR3M7ChqUo1gCqlwlNgAbQoWBeRjDHTgXTgV02tq4fwSrWwI0eOMGbMGNLS0hgzZgxH\njx496XqvvvoqaWlppKWl8eqrr/qe/+KLLzj77LPp27cvv//97323g/3tb39j0KBB2Gw2srPrnit8\n7LHH6Nu3L/369ePdd9/1Pd+mRngK7m1MBUBqrccpJ56rwxgzGpgPjBeRJu+F04tISrWwe++9l+jo\naDZv3sy2bdto164d27Zto1OnTr51jhw5Qnp6OnfffTdPPvkk+/btY+nSpdx2222cf/753HbbbSxZ\nsoQ9e/YwevRoVq9ezXfffUdJSQnjxo0jPj6efv368eabb7J//35GjRpF165dcblc7Nq1C4/Hw4ED\nB7jggguorKykY8eO7N69mzPPPJPt27e31K437yJS+3Qh3c+LSB82eRHJQfVFpFFUB86twHUi8k2t\ndYYCbwHjRGSnP9VqD1SpFrZmzRqKiooYMWIEaWlp7N27l65du3LJJZf4eqPvvvsuF198MUuWLOGu\nu+7C4XAwZ84cEhISyMvL47nnnuPuu++mXbt2rF27lnbt2vHcc8/x1ltv0b59exITE/noo4/o2rUr\nV1xxBbfccgtLly4lNzcXt9uNiDB48GC6deuGw+HgwQcfpKqqiu+++44ePXqQnt7qt1A2rSYX3p+l\nCSLiBm4H3gW+Bd4UkW+MMY8YY8afWO2/gHbA34wxXxlj1jZVrvZAlWphHTt2pFu3bqSlpfHBBx/4\nhkdMT0/nggsu4PHHH+eJJ55g8+bNfPHFFxQUFCAinHXWWRQWFuJ0OomLi6OsrAyPx0NiYiJ2u51u\n3bqRm5tLcXExdrudCy+8kO3btxMREUH//v3Zv38/OTk52O12rrvuOrp27cry5cspLy+nqqoKEWHC\nhAns3LmTVatWMXDgwGDvevN6oHHpwkA/e6DZobmRXnugSgXB6NGjGTx48M+WmjE8Dxw4wHvvvecb\nGrG8vJyvvvqK1atX+8r49ttvsdlsdO7cGWMMO3fupFu3btjtdsrKyoiNjUVEOHDgABEREbRv357D\nhw9js9nweDx88sknHDt2jLKyMjZv3ozD4cBut1NeXs6HH35I7969OXbsmG8MBI/Hw+rVq0lMTAzP\nsUZ1RHqlTg/vv/8+SUlJP3t+/vz5xMXF4Xa78Xq9HD58mMTERABKSkrIzc0FqkdfKiwsRESIiIgg\nPj4et9vNnj17qKysxOFw4HQ6fcE1JyeHPXv2ANUpwA6Hg/bt2+NyuSgtLcXpdDJ16lTfQDi7d+/m\n97//PQCTJk0iJSWFyMhIPB4PH330EX//+99b4V2yIMwHVNYAqlSQ1ATR3Nxcdu3a5VsOHjxIRUWF\n7+q50+mkffv2QPUoVf/xH//B2LFjKS0tpaioCJvNRmlpKXa73ddTrKioID4+niNHjpCQkIAxhr17\n9/py/auqqhg7diyAb3jCW265hYSEBGw2GzabjYqKClwuF59++ikVFRW+ofgiIyP54osvmDlzZmu/\nZY1rAyPS6zlQpVrY4cOH6d+/vy84Jicnc+jQIVwuFx6PB7vdzh133EFWVhYHDhzAGEP37t0pLCzE\nbrfXGZchMjLSN0vCsWPH6tRjs9mIjIzE6ay+qtKjRw8KCwuZOHEib7/9NhEREbjdbowxeL1eYmNj\nKS8vp2fPnuzbt4/ExEQ2bNjAsGHDgrXrzTsHGpUuJPt5DjRHz4EqdUrq0qULK1as8I1EdfjwYaKj\noxEREhIS6NSpE1lZWfTp04cuXbogIhQWFuJwOHy9yZqepsvlIiIigtLSUmJiYnwjdxljsNlsdOvW\njXbt2gEwZMgQUlJSWL16tW97m81GfHw8gG/0qH379hEbG8u0adOYPXt2670xTWkDPVANoEq1gtGj\nR9OlSxeioqJwOp0UFxcTERHBhRdeSJ8+fXA6neTn5xMTE+M7RK+5/ajmflG73U5ERAQVFRV4vV4S\nExN9vc6YmBjcbjf5+flERERgt9vZs2cPhw4dQkR8U9IYY7jyyiuJjIz03UIlIng8HjIyMiguLmb/\n/v2hfKvq0gCqlHI4HCxfvpyqqirsdju9evUiNjbWd0U9PT2dyy+/nISEBI4ePUpkZCSJiYlcfvnl\nvl5jzaH8wIEDSUhIoF+/fpx77rlUVVVRWVnJ0KFDGTRoEBMmTCAmJoaCgupEm4kTJ9KhQwfOPfdc\nzj77bI4ePUpVVRU9e/ZkyJAhdO3alc6dO3P//feTkpLi2y7k2sCAyhpAlWolV155JW+//TZer5eC\nggJ69eqF0+n0TZ3tcDjYvHkzZ599Nm63m5KSEjIyMrjkkkt8ty253W52797NhAkT2LVrFx07dmT4\n8OGICN9++y2TJk1i8+bNREVFUVpaSrdu3dizZw/Hjx/nyy+/ZOTIkYwbNw5jDIcPHyYnJ4cZM2aQ\nnJyM1+vF5QphNKpPb2NSStVWO4jWBLyaIFpSUuJL+awJom+88QYXX3yxL4jedtttuN1u/v73v/Pg\ngw/6gujixYtxu928+OKLviv6X3/9NWeeeSYej8d3YWrBggVMnjyZtLQ0Jk+ejNfrZdWqVcyZMwev\n18vBgwfDZ/K6NnAO1Hd+xM9FKRUEa9euFYfDIampqbJgwQI555xz5Le//a2sWbNGREQqKipkyJAh\n0r59exk+fLg8/fTTcvXVV4uIyB133CGRkZHSt29feemll6R3797idrtl3bp1kpaWJqmpqdK1a9c6\n9aWkpMiIESPk3HPPlaSkJGnfvr307dtXUlNTJT4+Xs4//3z57//+bxk+fHgwdzPQ+FJngfMEh/i3\nQHZz67PUxgA3UEoFSU3A69OnjyxatEhERB588ME6QTQzM1POPPNMGT58uOzevdu37aJFi6RPnz5y\n1llnyfr1633PT5kyRZKSksThcEhycrL89a9/FRGRoqIiGTlypPTt21dGjhwpGRkZsnbtWvF6vXLb\nbbdJnz59ZPDgwbJ169Zg7mLzA6gR/5YQBVC9D1Sp08xtt91GYmIiCxcubOmqmncfqEkX8PM+UEJz\nH6gOqKzUaeSVV15h7969LFu2LNRNOSVoAFXqNPHFF1/wxBNP8L//+7++NE7VPPouKnWaWLZsGUeO\nHOHSSy9lyJAhzJo1K9RNavP0HKhSqqU08xzoMIFP/Fw7Vs+BKqXUT2pSkcKXBlClVJgKbFrOUNAA\nqpQKU9oDVUopizSAKqWURQJUhLoRjdIAqpQKU3oOVCmlLNJDeKWUskh7oEopZZH2QJVSyiLtgSql\nlEVe9Cq8UkpZoofwSinVDHoIr5RSFmgPVCmlLNIAqpRSFulVeKWUskivwiullEV6CK+UUhaF/yG8\nTiqnlApTNT1Qf5amGWPGGWO+N8bsMsbMO8nrUcaYN068/pkxpldTZWoAVUqFqZoeqD9L44wxduDP\nwOXAQGCqMWZgvdVuBo6KSF/gKeDxpsrVAKqUClM1F5H8WZp0PrBLRPaISBWQBVxVb52rgFdP/P0W\nMMoY0+jMooGeA23WNKVKKeW//e/CggQ/V442xmTXerxcRJbXepwM5NV6nA9k1CvDt46IuI0xx4Au\nQFFDlepFJKVUWBKRcaFuQ1P0EF4pdTooAFJrPU458dxJ1zHGOIAOwOHGCtUAqpQ6HWwF0owxvY0x\nkcAUYG29ddYCM078nQl8ICLSWKF6CK+UOuWdOKd5O/AuYAdeEpFvjDGPANkishZ4EVhhjNkFHKE6\nyDbKNBFglVJKNUAP4ZVSyiINoEopZZEGUKWUskgDqFJKWaQBVCmlLNIAqpRSFmkAVUopi/4f3vsQ\nvZ1UNHwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f19ccb23080>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solving time step:  4\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVAAAADxCAYAAACd3+8mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl8VPW9+P/XZ5ashDWBhCSySGRV\nAYMRra2yCLaKLFFBqCiiLdIWr7Yu16Xoj6vY64al1ovXnZ9Ea6/AQxataO0tghLFWxEXlgSyABIg\nhCyTzPL+/hEyJpEkMyeTzIS8n4/HeTyYmfNZzpnhnc85n/P5fIyIoJRSKni2cFdAKaU6Kg2gSill\nkQZQpZSySAOoUkpZpAFUKaUs0gCqlFIWaQBVSp32jDEvGGO+M8bsaOJzY4x52hiz2xjzL2PM6EDy\n1QCqlOoMXgImN/P55UDGye0W4M+BZKoBVCl12hORfwBHm9nlKuAVqbUV6G6MSWkpX0ew9Qhyf6VU\n52Vak3iQMVIZ4L4H4EvAVe+tFSKyIojiUoGCeq8LT753oLlEwQZQpZRqF5XALwLcdzG4RCSzDatz\nShpAlVIRydCuAaoISK/3Ou3ke83Se6ARpkuXLuzdu/eUn7300kv86Ec/aucaKRUeNiA2wC0E1gLX\nn+yNvwA4LiLNXr7X1bHTee2118jMzKRLly6kpKRw+eWX889//tNyfsYYdu/e3eC9EydOcPvtt9O/\nf3/i4+M544wzyM7O5uOPP242r/LycgYOHGipHjU1NSxevJiMjAzi4+Pp378/8+bNIz8/31J+bWXx\n4sXMmTMn3NVoU/n5+Rhj8Hg84a5Kh2UAZ4Bbi3kZswrYAgw2xhQaY24yxvzSGPPLk7usB/YCu4Hn\ngFsDqWOnu4R/4oknWLp0Kc8++yyTJk0iKiqKjRs3smbNmqBbdx6PB4fjh6ewurqacePG0b17d95+\n+22GDh2Ky+Viw4YNbNiwgaysrIDzCkZ2djaFhYW89tprjBo1ioqKClauXMmmTZu46aabWpW3Cr1Q\nfOens1BewovIrBY+F2ChlYyD2Tq00tJSiY+PlzfeeKPJfT7++GO54IILpFu3bpKcnCwLFy6U6upq\n/+eALF++XAYNGiT9+/eXiy++WACJi4uT+Ph4ycnJkeeee06Sk5OlvLy82fo0zqvuvV27domISElJ\niVx55ZWSkJAgY8aMkfvuu08uuuiiU+b1t7/9TWJiYmT//v1NlldUVCRXXnml9OjRQ84880xZsWKF\n/7Pf//73kp2dLbNnz5YuXbrIiBEj5JtvvpGHH35YkpKSJC0tTd555x3//j/5yU/k7rvvljFjxkhC\nQoJMmTJFjhw5IiIiH3zwgaSmpjYou1+/fvK3v/1NNmzYIE6nUxwOh8THx8s555wjIrXfzbx58yQ5\nOVn69u0r9957r3g8nlMeR2VlpVx//fXSvXt3GTJkiDz66KMNyisqKpLp06dLYmKi9O/fX5YtW+b/\nzOVyyaJFiyQlJUVSUlJk0aJF4nK5GtT70UcflaSkJElOTpa33npL1q1bJxkZGdKjRw/5j//4D39e\nXq9XHnnkERk4cKD07NlTrr76av85SE9PF0Di4+MlPj5ePvroI3nxxRflwgsvlNtuu0169uwp99xz\nj/To0UP+9a9/+fM8dOiQxMbGynfffdfk99iBBBtfGmz9QJ4LcANyW1uela1TBdANGzaI3W4Xt9vd\n5D65ubmyZcsWcbvdkpeXJ0OGDJEnn3zS/zkgEyZMkCNHjkhlZaX/vbqgJyJy7bXXyty5c1usT0t5\nXXvttXL11VdLeXm5fPHFF9K3b98mA+hdd90lP/7xj5st7+KLL5YFCxZIVVWVbN++XRITE2XTpk0i\nUhtAo6OjZePGjeJ2u+XnP/+59O/fX5YsWSI1NTWyYsUKf5AXqQ2gffv2lS+++ELKy8tl+vTpMnv2\nbBFpPoDWlVW3b52pU6fKLbfcIuXl5XLo0CEZM2aMPPvss80e69GjR6WgoEDOPvtsf3ler1dGjx4t\nDz74oFRXV8uePXtkwIABsnHjRhERuf/++yUrK0sOHTok3333nYwdO1buu+8+f73tdrs8+OCD/mNO\nTEyUWbNmSVlZmezYsUNiYmJk7969IiLy1FNPSVZWlhQUFIjL5ZJbbrlFZs6cKSIieXl5AjT4rb34\n4otit9vl6aefFrfbLZWVlbJgwQK58847/fs89dRTcsUVVzT7PXYgrQpO/UFeDHDTANoOVq5cKX36\n9AkqzZNPPilTp071vwb8Qaf+e/UD6Pjx4+Wuu+7yv96+fbt069ZNEhIS5KyzzgooL4/HIw6HQ776\n6iv/Z/fcc0+TAXT+/Ply7bXXNnkc+/fvF5vNJmVlZf737r77bn+g//3vfy8TJkzwf7Z27VqJj4/3\ntwLLysoEkGPHjolIbQCtf4xffvmlOJ1O8Xg8QQfQgwcPSlRUlP+PiIjIa6+9Jpdccskpj6V+QBQR\nee655/zlbd26VdLT0xvs//DDD8sNN9wgIiIDBw6UdevW+T/buHGj9OvXT0RqA2hMTMwPjnnr1q3+\n/UePHi1vvfWWiIgMGTJE3nvvPf9nxcXF4nA4/H98TxVAG9etrr4+n09ERM477zx5/fXXT3ncHVCr\ngtMAkJUBbuEKoJ3qBkyvXr0oKSlp9t7Tt99+y+23305ubi6VlZV4PB7OO++8Bvukp6efMm39cg4c\n+L4Db+TIkZSWlvLee+8xf/78gPI6fPgwHo+nwef9+vVrtsxvv/22yc+Li4vp2bMnCQkJDfLLzc31\nv+7Tp4//37GxsSQmJmK32/2vobaTq3v37j+oe79+/XC73ZSUlDRZh6bs27cPt9tNSsr3Az98Pl+T\n56a4uLjBZ/X/vW/fPoqLi/11BPB6vVx88cX+tPXPY79+/SguLva/7tWr1w+OufF5KS8v95c1bdo0\nbLbv+2LtdjuHDh1q8lgbH1NWVhZxcXH8/e9/JyUlhd27dzNlypQm03cmdb3wkaxT9cKPHTuW6Oho\nVq9e3eQ+CxYsYMiQIezatYuysjIefvjh2qZ6PcY0P8Bi/PjxvPvuu1RUVLRYp6bySkpKwuFwUFDw\n/eCI/fv3N5nPhAkT+OSTTygsLDzl53379uXo0aOcOHGiQX6pqakt1rEpjevmdDpJTEwkPj6eysrv\nx5B4vV4OHz7sf934mNPT04mOjqakpITS0lJKS0spKyvjyy+/PGW5KSkpDY6zfj3S09MZMGCAP5/S\n0lJOnDjB+vXrgdrzsG/fvgb17tu3r6XjT09PZ8OGDQ3KcrlcpKamNvm9nur9uXPnsnLlSl599VWy\ns7OJiYmxVJ/TTSh74dtKpwqg3bp146GHHmLhwoWsXr2ayspK3G43GzZs4M477wRqHz/q2rUrXbp0\n4euvv+bPf255ToE+ffo0eHbz+uuvJyUlhWnTprFjxw68Xi8ul6tBa68ldrud6dOns3jxYiorK9m5\ncycvv/xyk/tPmDCBiRMnMm3aND799FM8Hg8nTpzg2Wef5YUXXiA9PZ0LL7yQe+65B5fLxb/+9S+e\nf/75Vj1OtHLlSnbu3EllZSUPPPAA2dnZ2O12zjrrLFwuF+vWrcPtdrNkyRKqq6v96fr06UN+fj4+\nnw+oDYiXXXYZd9xxB2VlZfh8Pvbs2cOHH354ynKvueYaHnnkEY4dO0ZRURHLly/3f3b++eeTkJDA\no48+SlVVFV6vlx07drBt2zYAZs2axZIlSzh8+DAlJSU89NBDls/BL3/5S+69915/QD58+DBr1qwB\nav8A2my2Jp/prW/OnDm89dZbrFy5kuuvv95SXU5XjgC3cOlUARTgjjvu4IknnmDJkiUkJSWRnp7O\n8uXLmTp1KgCPPfYYr732GgkJCdx8881ce+21Lea5ePFi5s6dS/fu3XnjjTeIiYnhgw8+YNiwYfzs\nZz+ja9euDB48mG3btvHGG28EXNfly5dTXl5OcnIyN9xwAzfeeGOz+7/55pv89Kc/5dprr6Vbt26M\nGDGC3NxcJkyYAMCqVavIz8+nb9++TJs2jQcffND/mRU///nPueGGG0hOTsblcvH0008DtX+onnnm\nGebPn09qairx8fGkpaX501199dVA7eXy6NG1s4a98sor1NTUMGzYMHr06EF2dnaD2yD1PfDAA6Sl\npTFgwAAmTJhAdnY20dHRQO0fnrfffpvPP/+cAQMGkJiYyPz58zl+/DgA9913H5mZmZxzzjmcffbZ\njB49mvvuu8/S8S9atIgpU6Zw2WWXkZCQwAUXXOB/zjcuLo57772Xiy66iO7du7N169Ym80lPT2f0\n6NEYY/y3GlTHaIGaxpenLdDJRBQAl1xyCXPmzPnBPd1w+POf/0xOTk6TLdaOYN68efTt25clS5aE\nuyqh1KrJRAYbI/8V4L6XwqeiY+GVatmBAwfYu3cvY8eOZdeuXTz++OP86le/Cne1LMvPz+d//ud/\n2L59e7irElG0E0mpNlBTU8MvfvELEhISGDduHFdddRW33hrQyLuIc//99zNixAh+97vfMWDAgHBX\nJ6LoJbxSqjNr1SX8cGNkVYD7nquX8Eop9b26Fmgk0wCqlIpI7TwfqCWRXj+lVCelLVCllLLIEPm9\n8BpAlVIRyQDOQCNUmOat1gCqlIpIxkDA801rAFVKqe8ZA057uGvRPH2QXoXEtm3bOOecc3C5XFRU\nVDB8+HB27NgR7mqpDqyuBRrIFrY66oP0KlTuu+8+XC4XVVVVpKWlcc8994S7Siq8WvUgfabTSG7P\nAAv6LjwP0msAVSFTU1PDmDFjiImJ4aOPPvJPTKw6rdYF0CgjuUkBFlSsI5FUB3fkyBHKy8txu924\nXC7i4+PDXSXVkXWAJ+m1BapCZsqUKcycOZO8vDwOHDjQYKJj1Sm1rgUabSQ3reX9AMxebYGqDuyV\nV17B6XRy3XXX4fV6ufDCC3n//fcZN25cuKumOioDRPhdIG2BKqXaSutaoLFGcgOc4c98pS1QpZT6\nngGiw12J5mkAVUpFpg7QiRTh1VNKdVoaQJVSqhUivBNJA6hSKjJpC1QppSzSAKqUUhZpL7xSSlnU\nAVqgOp1dK2zcuJHBgwczaNAgli5dGtK8CwoKuPTSSxk2bBjDhw9n2bJlIc2/Pq/Xy6hRo7jiiiva\nJP/S0lKys7MZMmQIQ4cOZcuWLSEv48knn2T48OGMGDGCWbNm4XK5Wp3nvHnz6N27NyNGjPC/d/To\nUSZOnEhGRgYTJ07k2LFjrS5HNaEugAayhYkGUIu8Xi8LFy5kw4YN7Ny5k1WrVrFz586Q5e9wOHj8\n8cfZuXMnW7du5U9/+lNI869v2bJlDB06tE3yBli0aBGTJ0/m66+/5v/+7/9CXlZRURFPP/00ubm5\n7NixA6/XS05OTqvzveGGG9i4cWOD95YuXcr48ePZtWsX48ePD/kfTlVP3VDOQLYw0QDaiNvtxuPx\n0NIQ108++YRBgwYxcOBAoqKimDlzJmvWrAlZPVJSUhg9ejQACQkJDB06lKKiopDlX6ewsJB169Yx\nf/78kOcNcPz4cf7xj39w0003ARAVFUX37t1DXo7H46GqqgqPx0NlZSV9+/ZtdZ4//vGP6dmz4YSU\na9asYe7cuQDMnTuX1atXt7oc1YQQt0CNMZONMd8YY3YbY+4+xednGGM+MMZsN8b8yxjz05by1ADa\niNfrJS8vr8UgWlRURHp6uv91WlpamwQ4gPz8fLZv305WVlbI877tttv4wx/+gM3WNj+FvLw8kpKS\nuPHGGxk1ahTz58+noqIipGWkpqby29/+ljPOOIOUlBS6devGZZddFtIy6hw6dIiUlBQAkpOTOXTo\nUJuUo/i+EymQraWsjLEDfwIuB4YBs4wxwxrtdh/whoiMAmYCz7SUrwbQUygsLCQ/P589e/a02BJt\na+Xl5cyYMYOnnnqKrl27hjTvt99+m969e3PeeeeFNN/6PB4Pn332GQsWLGD79u3Ex8eH/LL32LFj\nrFmzhry8PIqLi6moqGDlypUhLeNUjDEY06r5MlRzQtsCPR/YLSJ7RaQGyAGuarSPAHX/yboBxS1l\nqgG0CYWFhRQUFLB79258Pt8PPk9NTaWgoKDB/qmpqSGtg9vtZsaMGcyePZvp06eHNG+AzZs3s3bt\nWvr378/MmTN5//33mTNnTkjLSEtLIy0tzd96zs7O5rPPPgtpGe+99x4DBgwgKSkJp9PJ9OnT+eij\nj0JaRp0+ffpw4MABAA4cOEDv3r3bpBxFsAE00RiTW2+7pVFuqUBBvdeFJ9+rbzEwxxhTCKwHft1S\nFTWAtqCwsJA9e/bg8TRcN3XMmDHs2rWLvLw8ampqyMnJYcqUKSErV0S46aabGDp0KLfffnvI8q3v\nkUce8be2c3JyGDduXMhbbsnJyaSnp/PNN98AsGnTJoYNa3zl1DpnnHEGW7dupbKyEhFh06ZNbdYp\nNmXKFF5++WUAXn75Za66qnEjRoVU4AG0REQy620rLJQ2C3hJRNKAnwKvGmOajZEaQANQWFhIXl5e\ng0t6h8NBjx49mDRpEkOHDuWaa65h+PDhIStz8+bNvPrqqzz//POMHDmSkSNHsn79+pDlX9/kyZPb\nJN86f/zjHzn//PM555xz+Pzzz/n3f//3kOaflZVFdnY2o0ePpmvXrvh8Pm65pXEDJHizZs1i7Nix\nfPPNN6SlpfH8889z991385//+Z9kZGTw3nvvcffdP+iLUKES2l74IiC93uu0k+/VdxPwBoCIbAFi\ngMRmq6gTKjfkcrmafU7xyqnTqCg7Hlymdid43cGlcTjBE6Fp2rOs9kpj5TuykgZI6N6DsmNHg07X\nAbVuQuVkI7kB3lEyjzc/obIxxgF8C4ynNnBuA64TkS/r7bMBeF1EXjLGDAU2AanSTJCM8Of8I09F\n2XFYE+TfkasMrAsyzc8MbAoyzXgDW4NMc4GB/7Pwd/FcC+msprFyTFbOnZXvKNjfAnDiKu14CkgI\nh3KKiMcY8yvgHWrbrC+IyJfGmIeAXBFZC9wBPGeM+TdqG4s3NBc8QQOoUipShXgop4isp7ZzqP57\nD9T7907gomDy1ACqlIpMHWAsfIRXTynVaWkAVUqpVtAZ6ZVSygJtgSqllEU6oXLHc/jwYVwuF9HR\n0TrOWalw0hZox5OYmIjNZqOiooKYmBgcDj1FSoWFBtCOxxhDVFQUDocDl8uF2+0mJiZGW6NKtbe6\noZwRTANoE2w2G3FxcbjdbioqKoiOjsbpdNYO3wt2JIndUTtqJdg04y2kucBCmnMt/HGwks5qGivH\nZOXcWfmOrIwqsjuDT9MZaQu043M6nQ1ao3jd7Tcsc3OQaS6yOFRyzw+n62vRmTbYF+Q48H7O4Ms6\n02btmKycu/YY/gnBB+rOylA7nUcE09mYAmCMITY2lqioqHBXRZ0m2noGrNNCB1gTSVugQdAOJRUq\njRerU6fQAS7htQWqVITQZZRPQZc1VkoFQpdRbqQDXMJrAFUqQugyyo2EeFnjthDhdxiU6tw69TLK\nOpRTKRUqnW4ZZe1EUkq1RqdeRrkDXMJrAFUqgnXqZZQ7QADVVTkbaWlVznETJ1lYvdEBXk/L+3WU\nNAAOB3iCTGclTSSfB6vnzu5EPDU/eHvWrFn8/e9/p6SkhD59+vDggw8ydepUrrnmGvbv30+/fv14\n4403ftDRFMFatyrncCO5rwdY0NnNr8rZViL8DkME8rrbb1jml0GmGW6wHSwPKokvuQvdqg8EVw5w\nPDqFFNkbVJoDZmDQZR2PTrF0TFbOXbsM/4Qmx+mvWrXqlO9v2rQp+DJOBx3gHmiEV08p1WlpL7xS\nSlmkLVCllLJIA6hSSlmkAVQppawTnZFeKaWCJzaoifAJlTWANuJyufB6vdjtEf6nT6nTnBjw2AMd\n62NhVYUQ0ADaiMfjobq6Gp/Ph9PpxOl0YrPpgC2l2psYgzfgScx/ODChPWgAbaRLly7ExcUhIrjd\nbqqqqgD8wRS709qCZRdZSDM8yDQOR+1D5EGmOR6dElyak+kOmIFtX5aVY7Jy7qx+R8H+FgAcuqhc\noLwRfiWoAbQJdcsbR0VF4fP5cLvdVFZW1o5E2hrk6JMLrC32ZmUEjpXRQefLh0GlAfjE/ITL5a9B\npdlgZgRd1ifmJ5aOydLoJSuL1wX7W4DgVxntpASDN8LXNdZr0wDYbDaio6OJj48Pd1XUaUIXlWuZ\nYPBgD2gLFw2gSoVBIIvKPfnkkwwfPpwRI0Ywa9YsXC5XO9QscgiGGqID2sJFA6hSEaioqIinn36a\n3NxcduzYgdfrJScnJ9zVald1l/CBbOGiAVSpCOXxeKiqqsLj8VBZWUnfvn3DXaV2F8oAaoyZbIz5\nxhiz2xhzdxP7XGOM2WmM+dIY81pLeWonklIRKDU1ld/+9recccYZxMbGctlll3HZZZeFu1rtqu4e\naCgYY+zAn4CJQCGwzRizVkR21tsnA7gHuEhEjhljWpz+X1ugSkWgY8eOsWbNGvLy8iguLqaiooKV\nK1eGu1rtqvYS3hHQFoDzgd0isldEaoAcoPH0/jcDfxKRYwAi8l1LmWoAVSoCvffeewwYMICkpCSc\nTifTp0/no48+Cne12lVtJ1JUQBuQaIzJrbfd0ii7VKCg3uvCk+/VdxZwljFmszFmqzGmxUcl9BJe\nqQh0xhlnsHXrViorK4mNjWXTpk1kZrb7ihVhJRDMJXxJCJb0cAAZwCVAGvAPY8zZIlLaXAKlVITJ\nysoiOzub0aNH43A4GDVqFLfc0rhRdbozgV6eB6IISK/3Ou3ke/UVAh+LiBvIM8Z8S21A3dZkDXVR\nuYZaXFTuskngaYdF5dpr0TaHHTze4NIAxmFHgk1npSxLaSJ48ToAhxNxh2fsdjtr1ZCrIZnx8nzu\n0ID2/ZH5tNlF5YwxDuBbYDy1gXMbcJ2IfFlvn8nALBGZa4xJBLYDI0XkSFP5ags0WB63tSF/e4Kc\nLeZMm6UF2KwMlfy1/CGoNAB/NHeySqYGlWaWWR10WX80d1o6Jivnzsp3FPRvAWp/DyogoXrGU0Q8\nxphfAe8AduAFEfnSGPMQkCsia09+dpkxZifgBX7XXPAEDaBKqQgV6rHwIrIeWN/ovQfq/VuA209u\nAdEAqpSKSIKhOsKX5dQAqpSKSB1hNiYNoEqpiKQBVCmlWiGcU9UFQgOoUioiSWifA20TkV07pVSn\npZfwSillUW0vfFS4q9EsDaCNHD16lPLycowxp9yUUu2jI1zC61DORuqGcorIKbfJV1wZwUM5gx/2\naBw2xBP8mtp2h8HrCe7nYHPY8AVbVnsN5Wyv4Z8n002aMD6gZT06uFa1OPpn9pJ7c38W0L63mFeb\nHcrZViI7vIdRky1Oq0M59wUZdPs5La1GaWWlzGCHZELtsEy5Ibg05iWfpeGfVo7Jyrmz8h1ZHcoZ\nSPAsLS1l/vz57NixA2MML7zwAmPHjg2+vA5K74EqpSxbtGgRkydP5s0336SmpqZ2We1ORAOoUsqS\n48eP849//IOXXnoJgKioKKKiWu5QeeCBB+jZsye33XYbAPfeey+9e/dm0aJFbVndNtERhnLqjPRK\nRaC8vDySkpK48cYbGTVqFPPnz6eioqLFdPPmzeOVV14BwOfzkZOTw5w5c9q6um1CV+VUSlni8Xj4\n7LPPWLBgAdu3byc+Pp6lS5e2mK5///706tWL7du38+677zJq1Ch69erVDjVuG5EeQPUSXqkIlJaW\nRlpaGllZWQBkZ2cHFEAB5s+fz0svvcTBgweZN29eW1azTYVyVc62oi1QpSJQcnIy6enpfPPNNwBs\n2rSJYcOGBZR22rRpbNy4kW3btjFp0qS2rGabCvGqnG1CW6BKRag//vGPzJ49m5qaGgYOHMiLL74Y\nULqoqCguvfRSunfvjt0e2S24lmgvvFLKkpEjR5Kbmxt0Op/Px9atW/nLX/7SBrVqP3XLGkcyDaDB\ncjiDX9PG4ah96DrINAfMwKCSGIedDWZGUGlsDsMsszqoNAAOA+al4NLYLZRl5Zhw2IM+d5a+I7vD\n2vpGjiDLCcLOnTu54oormDZtGhkZGW1WTnvoCPdANYAGSxeVA3RROSAiF5UbNmwYe/cGNworUnWE\nsfCRXTulVKem90CVUsoCHcqplFIW6T1QpZSyqLYXPrLHwmsAVUpFJL2EV0qpVtAAqpRSFug9UKWU\nsqgjPAeqk4m0wOv1UlNTQ1VVVUDzMSoViMmTJ4e7ChGvbihnIFu46KJyjRw5coRPPvkEr9eLiGC3\n2/2bzWZj/KTJuqgctUNAfUEuKmeprNNxUTmHE3HXBLSr1+slMzOT1NRU3n777eDLCq9WDbnqljlI\nLsoNbOTaBjNDF5WLBDabDYfDQXR0dNOLym0N8u/IBcbS8E/bwfKgkviSu1haTC3YoZJQO1zSymJv\nVoZlWjkmK+fO0hDdYH8LUPt7CNCyZcsYOnQoZWVlwZdzGtBL+A4mNjYWh8Oha8CrsCssLGTdunXM\nnz8/3FUJi46wpEdkh3elOrHbbruNP/zhD5w4cSLcVQmLjvAcqLZAlYpAb7/9Nr179+a8884Ld1XC\nyoM9oC1ctAWqVATavHkza9euZf369bhcLsrKypgzZw4rV64Md9XajQ9bxA/l1BaoUhHokUceobCw\nkPz8fHJychg3blynCp51QnkP1Bgz2RjzjTFmtzHm7mb2m2GMEWNMi7362gJVSkWkUN4DNcbYgT8B\nE4FCYJsxZq2I7Gy0XwKwCPg4kHy1BapUhLvkkks64jOgrSaE9B7o+cBuEdkrIjVADnDVKfb7/4BH\nAVcgmWoAVUpFqKCWNU40xuTW225plFkqUFDvdeHJ974vzZjRQLqIrAu0hnoJr5SKSEFewpe0ZiSS\nMcYGPAHcEEw6DaDBcjiDGkkCWFu90eGoHR0TZJrgV6O084n5SXBpsL5aZtBlWVxhM+hzZ+U7sjuC\n/y1Am67KeToRDNWhG+deBKTXe5128r06CcAI4O8nB9EkA2uNMVNEpMm1pTWABsvjhk1BDt8bb2Bz\nkGkuMvBlkGmGWxv+GewKllC7iqWVIZZWVsu0NCzTwrmz9B0F+1uA2t+DalGIZ2PaBmQYYwZQGzhn\nAtf5yxI5DiTWvTbG/B34bXONCDjyAAARzUlEQVTBEzSAKqUiWKh64UXEY4z5FfAOYAdeEJEvjTEP\nAbkistZKvhpAlVIRKdRDOUVkPbC+0XsPNLHvJYHkqQFUKRWRBIPXF9lj4TWAKqUikvgM1a7IHsqp\nAVQpFZFEDF6PtkCVUip4ggbQ04nXG+TSEkopy0QMHndkB1AdyhkAj8dDRUUF1dXV4a6KOk20tKhc\nQUEBl156KcOGDWP48OEsW7asnWoWSQw+ryOgLVy0BdoMj8dDdXU1xhhiYmKw2yP7r6HqODZu3Njs\n5w6Hg8cff5zRo0dz4sQJzjvvPCZOnMiwYcPaqYYRQAC9hO9YRMQfOG02G7Gxsdhs9RrqdmfwI0ns\njtpRK8GmGd4+wz+PR6cEl+ZkOitDLIMuy+qwzGDPndXvyMqoogCGcqakpJCSUnuuEhISGDp0KEVF\nRZ0rgPoMuCI7REV27cKgpKQEt9v9w8BZx+uGdUEO3/uZhSF/Vod/WllZck/wyxpzpg32Bbm8cz9n\n8GWdabN2TO0xLHO8Cf63ALW/hyDk5+ezfft2srKygi+ro7OwanR70gDaSFJSErGxseGuhlIAlJeX\nM2PGDJ566im6du0a7uq0r9oJQSOaBlClIpTb7WbGjBnMnj2b6dOnh7s67U8DqFLKChHhpptuYujQ\nodx+++3hrk54CBDkXaL2po8xKRWBNm/ezKuvvsr777/PyJEjGTlyJOvXr2854elEgOoAtzDRFqhS\nEehHP/oRIhY6qE4negmvlFIWaQBVSimLNIAqpZRFGkBPQ3Zn0A9CWxqxYnVkjJWF0c600JfocNQ+\nGB9smmDLsnpM7TGqyO4I/rcAtb8hFRgNoKcZrxvWBHlz/yoLI1asjl7aGmSaCyyMXoLaoGZlhJCV\nNFaOqT1GFf3MBP9bgNrfg2qZD3CFuxLN0wCqlIpMegmvlFIWaQBVSimLNIAqpVQraABVSikLtAWq\nlFIW+YCqcFeieTqZSIBEhMrKynBXQ50mWloTCWqX/Rg8eDCDBg1i6dKl7VCrCCOAN8AtTDSABsDr\n9VJRUYHTqQ9Aq9BoaU0kr9fLwoUL2bBhAzt37mTVqlXs3LmznWoXQTwBbmGiAbQFbrebqqoq4uLi\nNICqkGmpBfrJJ58waNAgBg4cSFRUFDNnzmTNmjXtVLsIUXcPNIIDqN4DbYbL5cLn8xEfH48xBhEh\nLqErlcGOJLE0/NPC4nUOZ+0onGDTBDtU0mo6q2mCPSZLC/9Z/I6sjCpyOCksLGTy5MlNtkSLiopI\nT0/3v05LS+Pjjz8OvqyOTDuROiafz0dVVRUOh4O4uDig9h6ozWbj9df+f+x2O9HR0W1ah7ry27rV\nW1VVRVRUVJsv2VxVVUV0dPSpF+oLocrKSmJiYtq0nLr74VaXuv7d737Hvn37cLlczQbRTk+HcnY8\nx48f9//ncDhqT4/NZqOsrAyfz+dviXo8bfensW4iXZ/PR01NTZuVU1eGz2dhVc4giQhVVe3TpVpR\nUdHmgRpqg7WVchYvXszs2bMpKCigoKCAbt26MXbs2AaBNDU1lYKCAv/rwsJCUlNTQ1LvDkVboB1L\ndHQ0cXFx2Gw2RAS73e4PYnFxcW3eUvN6vVRXVxMbG4sxbTvphMfjwePxEBMT06blANTU1GCMaZf7\nyJWVlURHR7f5d1V3pWKlJbp69WruuusuSkpKcLlcfPfddw1ao2PGjGHXrl3k5eWRmppKTk4Or732\nWlscRuTqAJfw2onUSN3lX13wdLlcuFwuYmNj2/w/pIjgcrmIiYlp8+AJtQG0vTrG6lru7SE6Oprq\n6rZfKMdmsxEbG4vL5cLrDf5ZmkcffZSkpCTi4uIoKCigqKjI37nkcDhYvnw5kyZNYujQoVxzzTUM\nHz481IcQ2eoWlQtkCxMT5I/6tF+kpbKyki1btmC323G73VRXV7dbQHO73Rhj/LcO2lp1dXWb38ut\n4/P58Hq97Rawq6uriYqKapfvre5Wi9WW/N13382xY8coLS0lJSWF1NTU0+W+aKtOvknMFKbkBrbz\ni+ZTEclsTXlWaABtRET45z//aalFoZRVd955J6WlpURHRxMdHU1iYuLpEERbF0B7ZQo/CzCAvtpy\nADXGTAaWAXbgv0VkaaPPbwfmU3vj4DAwT0T2NZunBlClVBtpXQDtmSmMDzCAvtl8ADXG2IFvgYlA\nIbANmCUiO+vtcynwsYhUGmMWAJeIyLXNFav3QJVSkSm0QznPB3aLyF4RqQFygKsaFCfygYjUjdfe\nCqS1lKn2wiulIlNwvfCJxpj6zdUVIrKi3utUoKDe60Igq5n8bgI2tFSoBlClVGQKLoCWhKoTyRgz\nB8gEftLSvnoJr1QbO3r0KBMnTiQjI4OJEydy7NixU+738ssvk5GRQUZGBi+//LL//U8//ZSzzz6b\nQYMG8Zvf/Mb/ONhf/vIXhg8fjs1mIze34b3CRx55hEGDBjF48GDeeecd//sdaoan0D7GVASk13ud\ndvK9BowxE4B7gSki0uKzcNqJpFQbu/POO4mJiWHLli1s376dLl26sH37dnr06OHf5+jRo2RmZnLH\nHXfwxBNPsH//fpYtW8att97K+eefz6233srjjz/O3r17mTBhAqtXr+brr7+mrKyMyZMnk5CQwODB\ng3njjTc4cOAA48ePp3fv3rjdbnbv3o3X6+XgwYNccMEFVFdX0717d/bs2cOZZ57Jjh072urQW9eJ\n1DVTyAywE+mDFjuRHNR2Io2nNnBuA64TkS/r7TMKeBOYLCK7AilWW6BKtbE1a9ZQUlLC2LFjycjI\nYN++ffTu3ZtLLrnE3xp95513uPjii3n88ce5/fbbcTgcLFy4kMTERAoKCnjmmWe444476NKlC2vX\nrqVLly4888wzvPnmm3Tt2pWkpCQ+/PBDevfuzRVXXMHNN9/MsmXLyM/Px+PxICKMGDGCPn364HA4\nuP/++6mpqeHrr7+mb9++ZGa2+yOULasbCx/I1gIR8QC/At4BvgLeEJEvjTEPGWOmnNztP4EuwF+M\nMZ8bY9a2lK+2QJVqY927d6dPnz5kZGTw/vvv+6dHzMzM5IILLuDRRx/lscceY8uWLXz66acUFRUh\nIpx11lkUFxfjcrmIj4+noqICr9dLUlISdrudPn36kJ+fT2lpKXa7nQsvvJAdO3bgdDoZMmQIBw4c\nIC8vD7vdznXXXUfv3r1ZsWIFlZWV1NTUICJMnTqVXbt2sWrVKoYNGxbqQ29dCzQ+UxgWYAs0NzwP\n0msLVKkQmDBhAiNGjPjBVjeH58GDB3n33Xf9UyNWVlby+eefs3r1an8eX331FTabjZ49e2KMYdeu\nXfTp0we73U5FRQVxcXGICAcPHsTpdNK1a1eOHDmCzWbD6/WyefNmjh8/TkVFBVu2bMHhcGC326ms\nrOSDDz5gwIABHD9+3D8HgtfrZfXq1SQlJUXmXKM6I71SncN7771HcnLyD96/9957iY+Px+Px4PP5\nOHLkCElJSQCUlZWRn58P1M6+VFxcjIjgdDpJSEjA4/Gwd+9eqqurcTgcuFwuf3DNy8tj7969QO0Q\nYIfDQdeuXXG73ZSXl+NyuZg1a5Z/Ipw9e/bwm9/8BoDp06eTlpZGVFQUXq+XDz/8kL/+9a/tcJYs\niPAJlTWAKhUidUE0Pz+f3bt3+7dDhw5RVVXl7z13uVx07doVqJ2l6t/+7d+YNGkS5eXllJSUYLPZ\nKC8vx263+1uKVVVVJCQkcPToURITEzHGsG/fPv9Y/5qaGiZNmgTgn57w5ptvJjExEZvNhs1mo6qq\nCrfbzUcffURVVZV/Kr6oqCg+/fRT5s2b196nrHkdYEZ6vQeqVBs7cuQIQ4YM8QfH1NRUDh8+jNvt\nxuv1Yrfb+fWvf01OTg4HDx7EGENKSgrFxcXY7fYG8zJERUX5V0k4fvx4g3JsNhtRUVG4XLW9Kn37\n9qW4uJhp06bx1ltv4XQ68Xg8GGPw+XzExcVRWVlJv3792L9/P0lJSWzYsIHRo0eH6tBbdw80OlNI\nDfAeaJ7eA1XqtNSrVy9effVV/0xUR44cISYmBhEhMTGRHj16kJOTw8CBA+nVqxciQnFxMQ6Hw9+a\nrGtput1unE4n5eXlxMbG+mfuMsZgs9no06cPXbp0AWDkyJGkpaWxevVqf3qbzUZCQgKAf/ao/fv3\nExcXx+zZs1mwYEH7nZiWdIAWqAZQpdrBhAkT6NWrF9HR0bhcLkpLS3E6nVx44YUMHDgQl8tFYWEh\nsbGx/kv0useP6p4XtdvtOJ1Oqqqq8Pl8JCUl+VudsbGxeDweCgsLcTqd2O129u7dy+HDhxER/5I0\nxhiuvPJKoqKi/I9QiQher5esrCxKS0s5cOBAOE9VQxpAlVIOh4MVK1ZQU1OD3W6nf//+xMXF+XvU\nMzMzufzyy0lMTOTYsWNERUWRlJTE5Zdf7m811l3KDxs2jMTERAYPHsy5555LTU0N1dXVjBo1iuHD\nhzN16lRiY2MpKqodaDNt2jS6devGueeey9lnn82xY8eoqamhX79+jBw5kt69e9OzZ0/uuece0tLS\n/OnCrgNMqKwBVKl2cuWVV/LWW2/h8/koKiqif//+uFwu/9LZDoeDLVu2cPbZZ+PxeCgrKyMrK4tL\nLrnE/9iSx+Nhz549TJ06ld27d9O9e3fGjBmDiPDVV18xffp0tmzZQnR0NOXl5fTp04e9e/dy4sQJ\nPvvsM8aNG8fkyZMxxnDkyBHy8vKYO3cuqamp+Hw+3O4wRqPG9DEmpVR99YNoXcCrC6JlZWX+IZ91\nQfT111/n4osv9gfRW2+9FY/Hw1//+lfuv/9+fxBdunQpHo+H559/3t+j/8UXX3DmmWfi9Xr9HVOL\nFy9mxowZZGRkMGPGDHw+H6tWrWLhwoX4fD4OHToUOYvXdYB7oP77IwFuSqkQWLt2rTgcDklPT5fF\nixfLOeecI7/4xS9kzZo1IiJSVVUlI0eOlK5du8qYMWPkqaeekquvvlpERH79619LVFSUDBo0SF54\n4QUZMGCAeDweWbdunWRkZEh6err07t27QXlpaWkyduxYOffccyU5OVm6du0qgwYNkvT0dElISJDz\nzz9f/uu//kvGjBkTysMMNr402OA8wSGBbZDb2vIs1THIBEqpEKkLeAMHDpQlS5aIiMj999/fIIhm\nZ2fLmWeeKWPGjJE9e/b40y5ZskQGDhwoZ511lqxfv97//syZMyU5OVkcDoekpqbKf//3f4uISElJ\niYwbN04GDRok48aNk6ysLFm7dq34fD659dZbZeDAgTJixAjZtm1bKA+x9QHUSGBbmAKoPgeqVCdz\n6623kpSUxIMPPtjWRbXuOVCTKRDgc6CE5zlQnVBZqU7kpZdeYt++fSxfvjzcVTktaABVqpP49NNP\neeyxx/jf//1f/zBO1Tp6FpXqJJYvX87Ro0e59NJLGTlyJPPnzw93lTo8vQeqlGorrbwHOlpgc4B7\nx+k9UKWU+l7dUKTIpQFUKRWhgluWMxw0gCqlIpS2QJVSyiINoEopZZEAVeGuRLM0gCqlIpTeA1VK\nKYv0El4ppSzSFqhSSlmkLVCllLJIW6BKKWWRD+2FV0opS/QSXimlWkEv4ZVSygJtgSqllEUaQJVS\nyiLthVdKKYu0F14ppSzSS3illLIo8i/hdVE5pVSEqmuBBrK1zBgz2RjzjTFmtzHm7lN8Hm2Mef3k\n5x8bY/q3lKcGUKVUhKprgQayNc8YYwf+BFwODANmGWOGNdrtJuCYiAwCngQebSlfDaBKqQhV14kU\nyNai84HdIrJXRGqAHOCqRvtcBbx88t9vAuONMc2uLBrsPdBWLVOqlFKBO/AOLE4McOcYY0xuvdcr\nRGRFvdepQEG914VAVqM8/PuIiMcYcxzoBZQ0Vah2IimlIpKITA53HVqil/BKqc6gCEiv9zrt5Hun\n3McY4wC6AUeay1QDqFKqM9gGZBhjBhhjooCZwNpG+6wF5p78dzbwvohIc5nqJbxS6rR38p7mr4B3\nADvwgoh8aYx5CMgVkbXA88CrxpjdwFFqg2yzTAsBVimlVBP0El4ppSzSAKqUUhZpAFVKKYs0gCql\nlEUaQJVSyiINoEopZZEGUKWUsuj/AVlJsObcxLh7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f19cc941898>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solving time step:  5\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVAAAADxCAYAAACd3+8mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl4lOW5+PHvMzPZCWsCCUlkkciq\nAgYitlZkt0VkiQpCRQFtAVs82ioel6KHo9gjKpaqB4+Kyg8itUfgkkULWntUEKLYiriwJJAFkAAh\nZJlklvv3R8iYRJLMTCaZCdyf63qvi5l5n2UW7jzvcj+PERGUUkr5zhLsDiilVGulAVQppfykAVQp\npfykAVQppfykAVQppfykAVQppfykAVQpdd4zxrxijPneGLOnnteNMeY5Y8x+Y8y/jDGDvalXA6hS\n6kKwEhjXwOvXAalntzuBF7ypVAOoUuq8JyL/AE42sMsNwOtSZQfQ3hiT2Fi9Nl/74eP+SqkLl2lK\n4V7GSJmX+x6BrwB7jadWiMgKH5pLAnJrPM47+9yRhgr5GkCVUqpFlAG/8nLfRWAXkbRm7M45aQBV\nSoUkQ4sGqHwgpcbj5LPPNUjPgYaYNm3acPDgwXO+tnLlSn7605+2cI+UCg4LEOXlFgAbgFvPXo2/\nEjgtIg0evlf38YKzevVq0tLSaNOmDYmJiVx33XV89NFHftdnjGH//v21njtz5gz33HMP3bt3JyYm\nhosuuoiMjAw+/fTTBusqKSmhZ8+efvWjsrKSRYsWkZqaSkxMDN27d2fWrFnk5OT4VV9zWbRoETNm\nzAh2N5pVTk4OxhicTmewu9JqGSDMy63RuoxZA2wHehtj8owxs40xvzbG/PrsLpuAg8B+4CVgnjd9\nvOAO4Z9++mmWLFnCiy++yNixYwkPD2fLli2sX7/e59Gd0+nEZvvxR1hRUcGIESNo374977zzDn37\n9sVut7N582Y2b95Menq613X5IiMjg7y8PFavXs2gQYMoLS1l1apVbNu2jdmzZzepbhV4gfjOz2eB\nPIQXkWmNvC7AfH8q9mVr1YqKiiQmJkbWrl1b7z6ffvqpXHnlldKuXTtJSEiQ+fPnS0VFhed1QJYv\nXy69evWS7t27y9VXXy2AREdHS0xMjGRmZspLL70kCQkJUlJS0mB/6tZV/dy+fftERKSwsFCuv/56\niY2NlSFDhshDDz0kP/nJT85Z19/+9jeJjIyUw4cP19tefn6+XH/99dKhQwe5+OKLZcWKFZ7X/vCH\nP0hGRoZMnz5d2rRpIwMGDJBvv/1WHn/8cYmPj5fk5GR59913Pftfc801snDhQhkyZIjExsbKhAkT\n5MSJEyIi8sEHH0hSUlKttrt16yZ/+9vfZPPmzRIWFiY2m01iYmLksssuE5Gq72bWrFmSkJAgXbt2\nlQcffFCcTuc530dZWZnceuut0r59e+nTp488+eSTtdrLz8+XyZMnS1xcnHTv3l2WLVvmec1ut8uC\nBQskMTFREhMTZcGCBWK322v1+8knn5T4+HhJSEiQt99+WzZu3CipqanSoUMH+c///E9PXS6XS554\n4gnp2bOndOzYUW688UbPZ5CSkiKAxMTESExMjHzyySfy6quvylVXXSV33323dOzYUR544AHp0KGD\n/Otf//LUeezYMYmKipLvv/++3u+xFfE1vtTauoG85OUGZDW1PX+2CyqAbt68WaxWqzgcjnr3ycrK\nku3bt4vD4ZDs7Gzp06ePPPPMM57XARk1apScOHFCysrKPM9VBz0RkZtvvllmzpzZaH8aq+vmm2+W\nG2+8UUpKSuTLL7+Url271htA77//fvnZz37WYHtXX321zJ07V8rLy2X37t0SFxcn27ZtE5GqABoR\nESFbtmwRh8Mhv/zlL6V79+6yePFiqayslBUrVniCvEhVAO3atat8+eWXUlJSIpMnT5bp06eLSMMB\ntLqt6n2rTZw4Ue68804pKSmRY8eOyZAhQ+TFF19s8L2ePHlScnNz5dJLL/W053K5ZPDgwfLoo49K\nRUWFHDhwQHr06CFbtmwREZGHH35Y0tPT5dixY/L999/LsGHD5KGHHvL022q1yqOPPup5z3FxcTJt\n2jQpLi6WPXv2SGRkpBw8eFBERJ599llJT0+X3Nxcsdvtcuedd8rUqVNFRCQ7O1uAWr+1V199VaxW\nqzz33HPicDikrKxM5s6dK/fdd59nn2effVbGjx/f4PfYijQpOHUHedXLTQNoC1i1apV06dLFpzLP\nPPOMTJw40fMY8ASdms/VDKAjR46U+++/3/N49+7d0q5dO4mNjZVLLrnEq7qcTqfYbDb5+uuvPa89\n8MAD9QbQOXPmyM0331zv+zh8+LBYLBYpLi72PLdw4UJPoP/DH/4go0aN8ry2YcMGiYmJ8YwCi4uL\nBZBTp06JSFUArfkev/rqKwkLCxOn0+lzAD169KiEh4d7/oiIiKxevVqGDx9+zvdSMyCKiLz00kue\n9nbs2CEpKSm19n/88cfltttuExGRnj17ysaNGz2vbdmyRbp16yYiVQE0MjLyR+95x44dnv0HDx4s\nb7/9toiI9OnTR7Zu3ep5raCgQGw2m+eP77kCaN2+VffX7XaLiMgVV1whb7755jnfdyvUpODUA2SV\nl1uwAugFdQKmU6dOFBYWNnju6bvvvuOee+4hKyuLsrIynE4nV1xxRa19UlJSzlm2ZjtHjvxwAW/g\nwIEUFRWxdetW5syZ41Vdx48fx+l01nq9W7duDbb53Xff1ft6QUEBHTt2JDY2tlZ9WVlZnsddunTx\n/DsqKoq4uDisVqvnMVRd5Grfvv2P+t6tWzccDgeFhYX19qE+hw4dwuFwkJj4Q+KH2+2u97MpKCio\n9VrNfx86dIiCggJPHwFcLhdXX321p2zNz7Fbt24UFBR4Hnfq1OlH77nu51JSUuJpa9KkSVgsP1yL\ntVqtHDt2rN73Wvc9paenEx0dzd///ncSExPZv38/EyZMqLf8haT6Knwou6Cuwg8bNoyIiAjWrVtX\n7z5z586lT58+7Nu3j+LiYh5//PGqoXoNxjScYDFy5Ejee+89SktLG+1TfXXFx8djs9nIzf0hOeLw\n4cP11jNq1Ch27txJXl7eOV/v2rUrJ0+e5MyZM7XqS0pKarSP9anbt7CwMOLi4oiJiaGs7IccEpfL\nxfHjxz2P677nlJQUIiIiKCwspKioiKKiIoqLi/nqq6/O2W5iYmKt91mzHykpKfTo0cNTT1FREWfO\nnGHTpk1A1edw6NChWv3u2rWrX+8/JSWFzZs312rLbreTlJRU7/d6rudnzpzJqlWreOONN8jIyCAy\nMtKv/pxvAnkVvrlcUAG0Xbt2PPbYY8yfP59169ZRVlaGw+Fg8+bN3HfffUDV7Udt27alTZs2fPPN\nN7zwQuNzCnTp0qXWvZu33noriYmJTJo0iT179uByubDb7bVGe42xWq1MnjyZRYsWUVZWxt69e3nt\ntdfq3X/UqFGMHj2aSZMm8dlnn+F0Ojlz5gwvvvgir7zyCikpKVx11VU88MAD2O12/vWvf/Hyyy83\n6XaiVatWsXfvXsrKynjkkUfIyMjAarVyySWXYLfb2bhxIw6Hg8WLF1NRUeEp16VLF3JycnC73UBV\nQBwzZgz33nsvxcXFuN1uDhw4wIcffnjOdm+66SaeeOIJTp06RX5+PsuXL/e8NnToUGJjY3nyyScp\nLy/H5XKxZ88edu3aBcC0adNYvHgxx48fp7CwkMcee8zvz+DXv/41Dz74oCcgHz9+nPXr1wNVfwAt\nFku99/TWNGPGDN5++21WrVrFrbfe6ldfzlc2L7dguaACKMC9997L008/zeLFi4mPjyclJYXly5cz\nceJEAJ566ilWr15NbGwsd9xxBzfffHOjdS5atIiZM2fSvn171q5dS2RkJB988AH9+vXjF7/4BW3b\ntqV3797s2rWLtWvXet3X5cuXU1JSQkJCArfddhu33357g/u/9dZb/PznP+fmm2+mXbt2DBgwgKys\nLEaNGgXAmjVryMnJoWvXrkyaNIlHH33U85o/fvnLX3LbbbeRkJCA3W7nueeeA6r+UD3//PPMmTOH\npKQkYmJiSE5O9pS78cYbgarD5cGDq2YNe/3116msrKRfv3506NCBjIyMWqdBanrkkUdITk6mR48e\njBo1ioyMDCIiIoCqPzzvvPMOX3zxBT169CAuLo45c+Zw+vRpAB566CHS0tK47LLLuPTSSxk8eDAP\nPfSQX+9/wYIFTJgwgTFjxhAbG8uVV17puc83OjqaBx98kJ/85Ce0b9+eHTt21FtPSkoKgwcPxhjj\nOdWgWscI1NQ9PG2ETiaiABg+fDgzZsz40TndYHjhhRfIzMysd8TaGsyaNYuuXbuyePHiYHclkJo0\nmUhvY+S/vdz3WvhMNBdeqcYdOXKEgwcPMmzYMPbt28fSpUu56667gt0tv+Xk5PC///u/7N69O9hd\nCSl6EUmpZlBZWcmvfvUrYmNjGTFiBDfccAPz5nmVeRdyHn74YQYMGMDvf/97evToEezuhBQ9hFdK\nXciadAjf3xhZ4+W+l+shvFJK/aB6BBrKNIAqpUJSC88H6pdQ759S6gKlI1CllPKTIfSvwmsAVUqF\nJAOEeRuhgjRvtQZQpVRIMga8nm9aA6hSSv3AGAizBrsXDdMb6VVA7Nq1i8suuwy73U5paSn9+/dn\nz549we6WasWqR6DebEHro95IrwLloYcewm63U15eTnJyMg888ECwu6SCq0k30qeFGcnq6GVD3wfn\nRnoNoCpgKisrGTJkCJGRkXzyySeeiYnVBatpATTcSFa8lw0VaCaSauVOnDhBSUkJDocDu91OTExM\nsLukWrNWcCe9jkBVwEyYMIGpU6eSnZ3NkSNHak10rC5ITRuBRhjJSm58PwBzUEegqhV7/fXXCQsL\n45ZbbsHlcnHVVVfx/vvvM2LEiGB3TbVWBgjxs0A6AlVKNZemjUCjjGR5OcOf+VpHoEop9QMDRAS7\nEw3TAKqUCk2t4CJSiHdPKXXB0gCqlFJNEOIXkTSAKqVCk45AlVLKTxpAlVLKT3oVXiml/NQKRqA6\nnV0TbNmyhd69e9OrVy+WLFkS0Lpzc3O59tpr6devH/3792fZsmUBrb8ml8vFoEGDGD9+fLPUX1RU\nREZGBn369KFv375s37494G0888wz9O/fnwEDBjBt2jTsdnuT65w1axadO3dmwIABnudOnjzJ6NGj\nSU1NZfTo0Zw6darJ7ah6VAdQb7Yg0QDqJ5fLxfz589m8eTN79+5lzZo17N27N2D122w2li5dyt69\ne9mxYwd//vOfA1p/TcuWLaNv377NUjfAggULGDduHN988w3//Oc/A95Wfn4+zz33HFlZWezZsweX\ny0VmZmaT673tttvYsmVLreeWLFnCyJEj2bdvHyNHjgz4H05VQ3UqpzdbkGgArcPhcOB0OmksxXXn\nzp306tWLnj17Eh4eztSpU1m/fn3A+pGYmMjgwYMBiI2NpW/fvuTn5wes/mp5eXls3LiROXPmBLxu\ngNOnT/OPf/yD2bNnAxAeHk779u0D3o7T6aS8vByn00lZWRldu3Ztcp0/+9nP6Nix9oSU69evZ+bM\nmQDMnDmTdevWNbkdVY8Aj0CNMeOMMd8aY/YbYxae4/WLjDEfGGN2G2P+ZYz5eWN1agCtw+VykZ2d\n3WgQzc/PJyUlxfM4OTm5WQIcQE5ODrt37yY9PT3gdd9999388Y9/xGJpnp9CdnY28fHx3H777Qwa\nNIg5c+ZQWloa0DaSkpL43e9+x0UXXURiYiLt2rVjzJgxAW2j2rFjx0hMTAQgISGBY8eONUs7ih8u\nInmzNVaVMVbgz8B1QD9gmjGmX53dHgLWisggYCrwfGP1agA9h7y8PHJycjhw4ECjI9HmVlJSwpQp\nU3j22Wdp27ZtQOt+55136Ny5M1dccUVA663J6XTy+eefM3fuXHbv3k1MTEzAD3tPnTrF+vXryc7O\npqCggNLSUlatWhXQNs7FGIMxTZovQzUksCPQocB+ETkoIpVAJnBDnX0EqP5P1g4oaKxSDaD1yMvL\nIzc3l/379+N2u3/0elJSErm5ubX2T0pKCmgfHA4HU6ZMYfr06UyePDmgdQN8/PHHbNiwge7duzN1\n6lTef/99ZsyYEdA2kpOTSU5O9oyeMzIy+PzzzwPaxtatW+nRowfx8fGEhYUxefJkPvnkk4C2Ua1L\nly4cOXIEgCNHjtC5c+dmaUfhawCNM8Zk1djurFNbEpBb43He2edqWgTMMMbkAZuA3zTWRQ2gjcjL\ny+PAgQM4nbXXTR0yZAj79u0jOzubyspKMjMzmTBhQsDaFRFmz55N3759ueeeewJWb01PPPGEZ7Sd\nmZnJiBEjAj5yS0hIICUlhW+//RaAbdu20a9f3SOnprnooovYsWMHZWVliAjbtm1rtotiEyZM4LXX\nXgPgtdde44Yb6g5iVEB5H0ALRSStxrbCj9amAStFJBn4OfCGMabBGKkB1At5eXlkZ2fXOqS32Wx0\n6NCBsWPH0rdvX2666Sb69+8fsDY//vhj3njjDV5++WUGDhzIwIED2bRpU8Dqr2ncuHHNUm+1P/3p\nTwwdOpTLLruML774gn//938PaP3p6elkZGQwePBg2rZti9vt5s476w5AfDdt2jSGDRvGt99+S3Jy\nMi+//DILFy7kv/7rv0hNTWXr1q0sXPijaxEqUAJ7FT4fSKnxOPnsczXNBtYCiMh2IBKIa7CLOqFy\nbXa7vcH7FK+fNInS06d9qtOE2RCHs/EdW0mZlmzrfCsDENuhA8UnT/pcrhVq2oTKCUayvDyjZJY2\nPKGyMcYGfAeMpCpw7gJuEZGvauyzGXhTRFYaY/oC24AkaSBIhvh9/qGn9PRpusnXPpU5ZPpysfi2\nRvoBM4C+4tu5wq/NYAaKbzepf2GGMVQ+9KkMwE5zDcPkfZ/KbDcjfG5rp7nGr/fkz2fnz3fk628B\nqn4PygsBTOUUEacx5i7gXarGrK+IyFfGmMeALBHZANwLvGSM+TeqBou3NRQ8QQOoUipUBTiVU0Q2\nUXVxqOZzj9T4917gJ77UqQFUKRWaWkEufIh3Tyl1wdIAqpRSTaAz0iullB90BKqUUn7SCZVbn+PH\nj2O324mIiNA8Z6WCSUegrU9cXBwWi4XS0lIiIyOx2fQjUiooNIC2PsYYwsPDsdls2O12HA4HkZGR\nOhpVqqVVp3KGME3lrKNuKqfD4aCiooKIiAjCwsIYMXYM+Jq+Z7OC03X+lGnJts63MlSlgLorHT6X\na4WalsrZ00jWY1429MuGUzmbi45AGxEWFlZrNIrD6VfK3yXyT5/KfGcu9yuF0Z/0ypHyjk9lALaZ\n8YwR32bgf8/c4HNb28x4v96TP5+dP9+Rr78FqPo9KC8YqqbzCGE6G5MXjDFERUURHh4e7K6o80Rz\nz4B1XmgFayLpCNQHekFJBUrdxerUObSCi0g6AlUqROgyyuegyxorpbyhyyjX0QoO4TWAKhUidBnl\nOgK8rHFzCPEzDEpd2C7oZZQ1lVMpFSgX3DLKehFJKdUUF/Qyyq3gEF4DqFIh7IJeRrkVBFBN5ayj\nsVU5R44d4/tKjCGcWmhsVsSfdEQ/yrVUmVD+vKH+VM5p06bx97//ncLCQrp06cKjjz7KxIkTuemm\nmzh8+DDdunVj7dq1P7rQFMKalsrZ30jWm142dKmmcrYK4nC22GqZ/qQwXid/9anMZjOFW+Rln8oA\nrDazuU1e8KnMSjPX57ZWm9l+vaeWSv/09bcAVb+Hc1mzZs05n9+2bZvPbZwXWsE50BDvnlLqgqVX\n4ZVSyk86AlVKKT9pAFVKKT9pAFVKKf9JiM9IrwFUKRWSxAKVIT6hsgbQOux2Oy6XC6s1xP/0KXWe\nEwNOq7e5Pu5m7Ut9NIDW4XQ6qaiowO12ExYWRlhYGBaLJmwp1dLEGFxeT2Je2ax9qY9mItVRnYkk\nIjgcjqp1kMATTEeNGxuymUj+ZfpYEKfvf739KedfmfMsewldVM5bg9Is8n6Wd8fwHU25ZiKFkurl\njcPDw3G73TgcDsrKyhCHM6SzimbLcp/KvGzu4mH5d5/KAPyHeZxn5U6fytxtVvjc1n+Yx/16T6Ga\nvQRVvwfVOMHgCvF1jfXY1AsWi4WIiAhiYmKC3RV1ntBF5RonGJxYvdqCRQOoUkHgzaJyzzzzDP37\n92fAgAFMmzYNu93eAj0LHYKhkgivtmDRAKpUCMrPz+e5554jKyuLPXv24HK5yMzMDHa3WlT1Ibw3\nW7BoAFUqRDmdTsrLy3E6nZSVldG1a9dgd6nFBTKAGmPGGWO+NcbsN8YsrGefm4wxe40xXxljVjdW\np15EUioEJSUl8bvf/Y6LLrqIqKgoxowZw5gxY4LdrRZVfQ40EIwxVuDPwGggD9hljNkgIntr7JMK\nPAD8REROGWManf5fR6BKhaBTp06xfv16srOzKSgooLS0lFWrVgW7Wy2q6hDe5tXmhaHAfhE5KCKV\nQCZQd3r/O4A/i8gpABH5vrFKNYAqFYK2bt1Kjx49iI+PJywsjMmTJ/PJJ58Eu1stquoiUrhXGxBn\njMmqsdW9xy4JyK3xOO/sczVdAlxijPnYGLPDGNPorRJ6CK9UCLrooovYsWMHZWVlREVFsW3bNtLS\nWvw+8aAS8OUQvjAAN9LbgFRgOJAM/MMYc6mIFDVUQCkVYtLT08nIyGDw4MHYbDYGDRrEnXf6lrjQ\n+hlvD8+9kQ+k1HicfPa5mvKAT0XEAWQbY76jKqDuqreHmspZW6gsKtdSaZkWmwW3H6mcFpvB7fTt\n5+BPW+dd+ieayumtPmkx8nJWX6/2/an5rMFUTmOMDfgOGElV4NwF3CIiX9XYZxwwTURmGmPigN3A\nQBE5UV+9OgL1kTicDJUPfSqz01zDSHnHpzLbzHi/FmDzJ1XyHRnpUxmA8WYbJyXKpzIdTbnPbY03\n2/x6T/58dv58R77+FqDq96C8E6h7PEXEaYy5C3gXsAKviMhXxpjHgCwR2XD2tTHGmL2AC/h9Q8ET\nNIAqpUJUoHPhRWQTsKnOc4/U+LcA95zdvKIBVCkVkgRDRYgvy6kBVCkVklrDbEwaQJVSIUkDqFJK\nNUEwp6rzhgZQpVRIksDeB9osQrt3SqkLlh7CK6WUn6quwocHuxsN0gBax8mTJykpKcEYc85NKdUy\nWsMhvKZy1lFzVc5zbddd/wvEcf6kZVptBpePKZkANhs4fcxo9aetUE7/9DuV02ZlzMhRXi3r0co1\nacTRPa2TPJj1C6/2vdO8oatyhpL6RpzicPm1euMYWe9TmffMDdwmL/hUZqWZ69dKmb6mZEJVWmZl\nO9/KhJ8Wv9I//XlP/nx2/nxHvv4WoOr34E3wLCoqYs6cOezZswdjDK+88grDhl04K3rqOVCllN8W\nLFjAuHHjeOutt6isrKSsrCzYXWpRGkCVUn45ffo0//jHP1i5ciUA4eHhhIc3fkHlkUceoWPHjtx9\n990APPjgg3Tu3JkFCxY0Z3ebRWtI5dQZ6ZUKQdnZ2cTHx3P77bczaNAg5syZQ2lpaaPlZs2axeuv\nvw6A2+0mMzOTGTNmNHd3m4WuyqmU8ovT6eTzzz9n7ty57N69m5iYGJYsWdJoue7du9OpUyd2797N\ne++9x6BBg+jUqVML9Lh5hHoA1UN4pUJQcnIyycnJpKenA5CRkeFVAAWYM2cOK1eu5OjRo8yaNas5\nu9msArkqZ3PREahSISghIYGUlBS+/fZbALZt20a/fv28Kjtp0iS2bNnCrl27GDt2bHN2s1kFeFXO\nZqEjUKVC1J/+9CemT59OZWUlPXv25NVXX/WqXHh4ONdeey3t27fHag3tEVxj9Cq8UsovAwcOJCsr\ny+dybrebHTt28Je//KUZetVyqpc1DmWaiVRH44vKjQ7hTCTfF3rzJ6MIqv7y+lrMn7b8eU8hvRAd\nzbuo3N69exk/fjyTJk1i6dKlzdKGD5qUiZSQliQzsuZ6te9S87BmIrUG4nDponLoonIQmovK9evX\nj4MHDzZb/S2pNeTCh3bvlFIXND0HqpRSftBUTqWU8lNruA9UA6hSKiRVXYUP7Vx4DaBKqZCkh/BK\nKdUEGkCVUsoPeg5UKaX81BruA9XJRBrhcrmorKykvLzcq/kYlfLGuHHjgt2FkFedyunNFiyaylnH\niRMn2LlzJy6XCxHBarV6NovFwqhxY0I4lVMXlYNWsKicD6mcLpeLtLQ0kpKSeOcd3zKlQkCTUjnb\npfWSn2T90at9N5spmsoZCiwWCzabjYiIiHoXlRso9efKn8sXZphfC9FdJ3/1qcxmM4XZstynMi+b\nu3xOlYSqdEl/FnvzJy3Tn/fkz2fnz3fk628Bqn4P3lq2bBl9+/aluLjY53bOB3oI38pERUVhs9l0\nDXgVdHl5eWzcuJE5c+YEuytB0RqW9Ajt8K7UBezuu+/mj3/8I2fOnAl2V4KiNdwHqiNQpULQO++8\nQ+fOnbniiiuC3ZWgcmL1agsWHYEqFYI+/vhjNmzYwKZNm7Db7RQXFzNjxgxWrVoV7K61GDeWkE/l\n1BGoUiHoiSeeIC8vj5ycHDIzMxkxYsQFFTyrBfIcqDFmnDHmW2PMfmPMwgb2m2KMEWNMo1f1dQSq\nlApJgTwHaoyxAn8GRgN5wC5jzAYR2Vtnv1hgAfCpN/XqCFSpEDd8+PDWeA9okwkBPQc6FNgvIgdF\npBLIBG44x37/ATwJ2L2pVAOoUipE+bSscZwxJqvGVvcm5SQgt8bjvLPP/dCaMYOBFBHZ6G0P9RBe\nKRWSfDyEL2xKJpIxxgI8DdzmUzlN5ayt8VU5xyAOH3MY/Uj5C+X0z6pyvq+WGcppmaGaytnKNSkb\nJSJtgCRmveXVvodM3wZTOY0xw4BFIjL27OMHAETkibOP2wEHgJKzRRKAk8AEEal3bWkdgfpIHE76\nyuc+lfnaDA7p9E9fV7CEqlUsb5MXfCqz0sz1a7XMUE3L/MIM8/m3AFW/B9W4AM/GtAtINcb0APKB\nqcAtnrZETgNx1Y+NMX8HftdQ8AQNoEqpEBaoq/Ai4jTG3AW8C1iBV0TkK2PMY0CWiGzwp14NoEqp\nkBToVE4R2QRsqvPcI/XsO9yq15CFAAARjUlEQVSbOjWAKqVCkmBwuUM7F14DqFIqJInbUGEP7VRO\nDaBKqZAkYnA5dQSqlFK+EzSAnk9cLt/v+VNK+UfE4HSEdgDVVE4vOJ1OSktLqaioCHZX1HmisUXl\ncnNzufbaa+nXrx/9+/dn2bJlLdSzUGJwu2xebcGiI9AGOJ1OKioqMMYQGRmJ1Rrafw1V67Fly5YG\nX7fZbCxdupTBgwdz5swZrrjiCkaPHk2/fv1aqIchQIAQP4TXVM46ysvL+eijj6ioqMBisRAREYHF\n8sNAvaVSOUM5/dPfcuddWmYLpnLecMMN3HXXXYwePdrn9oKoSamcZkCa8L8NJgL9oLfRVTlDQWFh\nIQ6Hg6ioqFqBs5o4nFwse3yq84AZwCXyT5/KfGcub7H0z5Hi+1Rp28x4xsh6n8q8Z27wua1tZnyL\npWX68x35+luAqt+DL3Jycti9ezfp6ek+t9Xq+ThWaWkaQOuIj48nKioq2N1QCoCSkhKmTJnCs88+\nS9u2bYPdnZZVNSFoSNMAqlSIcjgcTJkyhenTpzN58uRgd6flaQBVSvlDRJg9ezZ9+/blnnvuCXZ3\ngkOAEJ/1T29jUioEffzxx7zxxhu8//77DBw4kIEDB7Jp06bGC55PBKjwcgsSHYEqFYJ++tOf4uMd\nMucfPYRXSik/aQBVSik/aQBVSik/tYIAqplIdYTKonKhnL3kb7nzLqtIF5VrTNMykXqlCX/0MhNp\nimYitQricNJNvvapzCHT16/spZZavG6ofOhTGYCd5hq/MoR8bWunuaZFFnv72gz26zvy9bcAVb8H\n5QU3YA92JxqmAVQpFZpawSG8BlClVGjSAKqUUn7SAKqUUk2gAVQppfygI1CllPKTGygPdicappOJ\neElEKCsrC3Y31HmisTWRoGrZj969e9OrVy+WLFnSAr0KMQK4vNyCRAOoF1wuF6WlpYSFhQW7K+o8\n0diaSC6Xi/nz57N582b27t3LmjVr2Lt3bwv1LoQ4vdyCRANoIxwOB+Xl5URHR2sAVQHT2Ah0586d\n9OrVi549exIeHs7UqVNZv963JVRavepzoCEcQPUcaAPsdjtut5uYmBiMMYgI0W1jfc4kMWE2n9fB\nMWE2vjaDfS7zhRnmc5md5hqfylSVs7LdjGj2tvx9T/58dv58R/5kFZkwG3l5eYwbN67ekWh+fj4p\nKSmex8nJyXz66ac+t9Wq6UWk1sntdlNeXo7NZiM6OhqoOgdqsVh48/+txmq1EhER0ax9qG6/uUe9\n5eXlhIeHN/uSzeXl5T9a4bQ5lJWVERkZ2aztVJ8P93ep69///vccOnQIu93eYBC94GkqZ+tz+vRp\nz38Om63q47FYLBQXF+N2uz0jUaez+f40Vk/w4na7qaysbLZ2qttwu31f1thXIkJ5ectcUi0tLW32\nQA1VwdqfdhYtWsT06dPJzc0lNzeXdu3aMWzYsFqBNCkpidzcXM/jvLw8kpKSAtLvVkVHoK1LREQE\n0dHRWCwWRASr1eoJYtHR0c0+UnO5XFRUVBAVFYUxTZrMplFOpxOn00lkZGSztgNQWVmJMaZFziOX\nlZURERHR7N9V9ZGKPyPRdevWcf/991NYWIjdbuf777+vNRodMmQI+/btIzs7m6SkJDIzM1m9enVz\nvI3Q1QoO4fUiUh3Vh3/VwdNut2O324mKimr2/5Aigt1uJzIystmDJ1QF0Ja6MFY9cm8JERERVFQ0\n/0I5FouFqKgo7HY7Lpfv99I8+eSTxMfHEx0dTW5uLvn5+Z6LSzabjeXLlzN27Fj69u3LTTfdRP/+\n/QP9FkJb9aJy3mxBovOB1lFWVsb27duxWq04HA4qKipaLKA5HA6MMZ5TB82toqKi2c/lVnO73bhc\nrhYL2BUVFYSHh7fI91Z9qsXfkfzChQs5deoURUVFJCYmkpSUdL6cF23afKBxacIEL+cDfTU484Fq\nAK1DRPjoo4/8GlEo5a/77ruPoqIiIiIiiIiIIC4u7nwIok0LoJ3ShF94GUDfaDyAGmPGAcsAK/A/\nIrKkzuv3AHOoOnFwHJglIocarFMDqFKqmTQtgHZME0Z6GUDfajiAGmOswHfAaCAP2AVME5G9Nfa5\nFvhURMqMMXOB4SJyc0PN6jlQpVRoCmwq51Bgv4gcFJFKIBO4oVZzIh+ISHW+9g4gubFK9Sq8Uio0\n+XYVPs4YU3O4ukJEVtR4nATk1nicB6Q3UN9sYHNjjWoAVUqFJt8CaGGgLiIZY2YAaUCjaXN6CK9U\nMzt58iSjR48mNTWV0aNHc+rUqXPu99prr5Gamkpqaiqvvfaa5/nPPvuMSy+9lF69evHb3/7WczvY\nX/7yF/r374/FYiErq/a5wieeeIJevXrRu3dv3n33Xc/zrWqGp8DexpQPpNR4nHz2uVqMMaOAB4EJ\nItLovXB6EUmpZnbfffcRGRnJ9u3b2b17N23atGH37t106NDBs8/JkydJS0vj3nvv5emnn+bw4cMs\nW7aMefPmMXToUObNm8fSpUs5ePAgo0aNYt26dXzzzTcUFxczbtw4YmNj6d27N2vXruXIkSOMHDmS\nzp0743A42L9/Py6Xi6NHj3LllVdSUVFB+/btOXDgABdffDF79vi2GqkPmnYRqW2akOblRaQPGr2I\nZKPqItJIqgLnLuAWEfmqxj6DgLeAcSKyz5tmdQSqVDNbv349hYWFDBs2jNTUVA4dOkTnzp0ZPny4\nZzT67rvvcvXVV7N06VLuuecebDYb8+fPJy4ujtzcXJ5//nnuvfde2rRpw4YNG2jTpg3PP/88b731\nFm3btiU+Pp4PP/yQzp07M378eO644w6WLVtGTk4OTqcTEWHAgAF06dIFm83Gww8/TGVlJd988w1d\nu3YlLa3Fb6FsXHUuvDdbI0TECdwFvAt8DawVka+MMY8ZYyac3e2/gDbAX4wxXxhjNjRWr45AlWpm\n7du3p0uXLqSmpvL+++97pkdMS0vjyiuv5Mknn+Spp55i+/btfPbZZ+Tn5yMiXHLJJRQUFGC324mJ\niaG0tBSXy0V8fDxWq5UuXbqQk5NDUVERVquVq666ij179hAWFkafPn04cuQI2dnZWK1WbrnlFjp3\n7syKFSsoKyujsrISEWHixIns27ePNWvW0K9fv0C/9aaNQGPShH5ejkCzgnMjvY5AlQqAUaNGMWDA\ngB9t1XN4Hj16lPfee88zNWJZWRlffPEF69at89Tx9ddfY7FY6NixI8YY9u3bR5cuXbBarZSWlhId\nHY2IcPToUcLCwmjbti0nTpzAYrHgcrn4+OOPOX36NKWlpWzfvh2bzYbVaqWsrIwPPviAHj16cPr0\nac8cCC6Xi3Xr1hEfHx+ac43qjPRKXRi2bt1KQkLCj55/8MEHiYmJwel04na7OXHiBPHx8QAUFxeT\nk5MDVM2+VFBQgIgQFhZGbGwsTqeTgwcPUlFRgc1mw263e4JrdnY2Bw8eBKpSgG02G23btsXhcFBS\nUoLdbmfatGmeiXAOHDjAb3/7WwAmT55McnIy4eHhuFwuPvzwQ/7617+2wKfkhxCfUFkDqFIBUh1E\nc3Jy2L9/v2c7duwY5eXlnqvndrudtm3bAlWzVP3bv/0bY8eOpaSkhMLCQiwWCyUlJVitVs9Isby8\nnNjYWE6ePElcXBzGGA4dOuTJ9a+srGTs2LEAnukJ77jjDuLi4rBYLFgsFsrLy3E4HHzyySeUl5d7\npuILDw/ns88+Y9asWS39kTWsFcxIr+dAlWpmJ06coE+fPp7gmJSUxPHjx3E4HLhcLqxWK7/5zW/I\nzMzk6NGjGGNITEykoKAAq9Vaa16G8PBwzyoJp0+frtWOxWIhPDwcu73qqkrXrl0pKChg0qRJvP32\n24SFheF0OjHG4Ha7iY6OpqysjG7dunH48GHi4+PZvHkzgwf7Npt/A5p2DjQiTUjy8hxotp4DVeq8\n1KlTJ9544w3PTFQnTpwgMjISESEuLo4OHTqQmZlJz5496dSpEyJCQUEBNpvNM5qsHmk6HA7CwsIo\nKSkhKirKM3OXMQaLxUKXLl1o06YNAAMHDiQ5OZl169Z5ylssFmJjYwE8s0cdPnyY6Ohopk+fzty5\nc1vug2lMKxiBagBVqgWMGjWKTp06ERERgd1up6ioiLCwMK666ip69uyJ3W4nLy+PqKgozyF69e1H\n1feLWq1WwsLCKC8vx+12Ex8f7xl1RkVF4XQ6ycvLIywsDKvVysGDBzl+/Dgi4lmSxhjD9ddfT3h4\nuOcWKhHB5XKRnp5OUVERR44cCeZHVZsGUKWUzWZjxYoVVFZWYrVa6d69O9HR0Z4r6mlpaVx33XXE\nxcVx6tQpwsPDiY+P57rrrvOMGqsP5fv160dcXBy9e/fm8ssvp7KykoqKCgYNGkT//v2ZOHEiUVFR\n5OdXJdpMmjSJdu3acfnll3PppZdy6tQpKisr6datGwMHDqRz58507NiRBx54gOTkZE+5oGsFEypr\nAFWqhVx//fW8/fbbuN1u8vPz6d69O3a73bN0ts1mY/v27Vx66aU4nU6Ki4tJT09n+PDhntuWnE4n\nBw4cYOLEiezfv5/27dszZMgQRISvv/6ayZMns337diIiIigpKaFLly4cPHiQM2fO8PnnnzNixAjG\njRuHMYYTJ06QnZ3NzJkzSUpKwu1243AEMRrVpbcxKaVqqhlEqwNedRAtLi72pHxWB9E333yTq6++\n2hNE582bh9Pp5K9//SsPP/ywJ4guWbIEp9PJyy+/7Lmi/+WXX3LxxRfjcrk8F6YWLVrElClTSE1N\nZcqUKbjdbtasWcP8+fNxu90cO3YsdBavawXnQD3nR7zclFIBsGHDBrHZbJKSkiKLFi2Syy67TH71\nq1/J+vXrRUSkvLxcBg4cKG3btpUhQ4bIs88+KzfeeKOIiPzmN7+R8PBw6dWrl7zyyivSo0cPcTqd\nsnHjRklNTZWUlBTp3LlzrfaSk5Nl2LBhcvnll0tCQoK0bdtWevXqJSkpKRIbGytDhw6V//7v/5Yh\nQ4YE8m36Gl9qbXCFYBPvNshqant+9dHHAkqpAKkOeD179pTFixeLiMjDDz9cK4hmZGTIxRdfLEOG\nDJEDBw54yi5evFh69uwpl1xyiWzatMnz/NSpUyUhIUFsNpskJSXJ//zP/4iISGFhoYwYMUJ69eol\nI0aMkPT0dNmwYYO43W6ZN2+e9OzZUwYMGCC7du0K5FtsegA14t0WpACq94EqdYGZN28e8fHxPPro\no83dVNPuAzVpAl7eB0pw7gPVCZWVuoCsXLmSQ4cOsXz58mB35bygAVSpC8Rnn33GU089xf/93/95\n0jhV0+inqNQFYvny5Zw8eZJrr72WgQMHMmfOnGB3qdXTc6BKqebSxHOggwU+9nLvaD0HqpRSP6hO\nRQpdGkCVUiHKt2U5g0EDqFIqROkIVCml/KQBVCml/CRAebA70SANoEqpEKXnQJVSyk96CK+UUn7S\nEahSSvlJR6BKKeUnHYEqpZSf3OhVeKWU8osewiulVBPoIbxSSvlBR6BKKeUnDaBKKeUnvQqvlFJ+\n0qvwSinlJz2EV0opP4X+IbwuKqeUClHVI1BvtsYZY8YZY741xuw3xiw8x+sRxpg3z77+qTGme2N1\nagBVSoWo6hGoN1vDjDFW4M/AdUA/YJoxpl+d3WYDp0SkF/AM8GRj9WoAVUqFqOqLSN5sjRoK7BeR\ngyJSCWQCN9TZ5wbgtbP/fgsYaYxpcGVRX8+BNmmZUqWU8t6Rd2FRnJc7Rxpjsmo8XiEiK2o8TgJy\nazzOA9Lr1OHZR0ScxpjTQCegsL5G9SKSUiokici4YPehMXoIr5S6EOQDKTUeJ5997pz7GGNsQDvg\nREOVagBVSl0IdgGpxpgexphwYCqwoc4+G4CZZ/+dAbwvItJQpXoIr5Q67509p3kX8C5gBV4Rka+M\nMY8BWSKyAXgZeMMYsx84SVWQbZBpJMAqpZSqhx7CK6WUnzSAKqWUnzSAKqWUnzSAKqWUnzSAKqWU\nnzSAKqWUnzSAKqWUn/4/zhElqK3sXRwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f19ccb13400>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# define iteration parameters\n",
    "newton_tol = 1e-6\n",
    "t = .0\n",
    "T = 1\n",
    "k = 0\n",
    "times = [t]\n",
    "\n",
    "# Time loop\n",
    "while t < T:\n",
    "    # Increment time\n",
    "    t += dt\n",
    "    k += 1\n",
    "    times.append(t)\n",
    "    p0 = p.val\n",
    "    print('Solving time step: ', k)\n",
    "    # solve newton iteration\n",
    "    err = np.inf\n",
    "    while err > newton_tol:\n",
    "        eq = f(p, p0)   \n",
    "        p = p - sps.linalg.spsolve(eq.jac, eq.val)\n",
    "        err = np.sqrt(np.sum(eq.val**2))\n",
    "    pp.plot_grid(g, p.val,color_map = [0, 1])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
